Chiller performance test working condition self-optimization method and system based on cross-model knowledge embedding

Abstract

The application discloses a chiller performance test working condition self-optimization method and system based on cross-model knowledge embedding. The method calls a cross-model knowledge embedding mapping structure pre-trained through meta-task learning, fuses inherent attribute information of a to-be-tested unit and initial working conditions across domains by using a unit structure coding path, a working condition condition coding path and a self-adaptive kernel form description quantity fusion path, generates an initial kernel form description quantity set, constructs a Gaussian process proxy model based on the set and generates a performance proxy surface, combines a refrigeration cycle thermodynamic constraint boundary, performs extreme value search on a joint evaluation function of fusion information collection benefits and feasible domain risk penalties, determines a next test working condition, drives the unit to perform testing and incrementally updates the proxy surface, and iterates repeatedly until global prediction uncertainty meets a termination condition. The method can quickly and accurately depict global performance of the chiller with a small number of test working conditions, and significantly improves test efficiency and reduces test cost.

Description

Technical Field

[0001] This invention relates to the field of data processing, and in particular to a method and system for self-optimization of chiller performance testing conditions based on cross-model knowledge embedding. Background Technology

[0002] Chiller performance testing is a crucial step in evaluating the energy efficiency of chillers under different operating conditions. Its purpose is to obtain a comprehensive picture of the unit's performance response across the entire operating range, providing a foundation for unit selection, matching, and system optimization. Current chiller performance testing typically relies on a pre-set standard operating condition table, stably operating the chiller under several discrete combinations of cooling water temperature and load rate, recording the corresponding performance, and then obtaining the unit's performance curve through data fitting. This fixed-condition traversal testing method organizes the test sequence independently of the actual performance characteristics of the unit gradually revealed during the test. It is difficult to dynamically adjust the subsequent test direction based on the acquired performance information, resulting in a significant waste of testing resources on redundant operating conditions that contribute little to improving the performance model, while the key areas that truly determine the performance curve shape may be undersampled because they are not included in the pre-set operating condition table. Therefore, how to adaptively plan the test condition sequence to achieve an accurate characterization of the chiller's performance across all operating conditions within limited test time and energy consumption has become an urgent problem to be solved. Summary of the Invention

[0003] This invention provides a self-optimization method and system for chiller performance testing based on cross-model knowledge embedding.

[0004] In a first aspect, embodiments of the present invention provide a self-optimization method for chiller performance testing conditions based on cross-model knowledge embedding, including: The pre-built cross-aircraft knowledge embedding mapping structure is retrieved. The cross-aircraft knowledge embedding mapping structure is obtained by learning and training through meta-tasks from historical multi-aircraft test sequences. It includes the operating condition encoding path for mapping the operating condition sequence to the operating condition embedding representation, the unit structure encoding path for mapping the unit's inherent attribute information to the structure embedding representation, and the adaptive kernel morphology descriptor fusion path for fusing the operating condition embedding representation and the structure embedding representation and outputting the kernel morphology descriptor set. The initial performance observations of the chiller under test are obtained under the initial test conditions. The inherent attribute information of the chiller under test is retrieved. The inherent attribute information is used as input to generate a structural embedding representation through the unit structure coding path. The initial test conditions are used as input to generate a condition embedding representation through the condition coding path. The structural embedding representation and the condition embedding representation are fused across domains through an adaptive kernel morphological descriptor fusion path. The initial kernel morphological descriptor set is output, and the initial performance observations are used as the initial observation records. A Gaussian process surrogate model is constructed based on the initial set of kernel morphology descriptors. The posterior distribution of the Gaussian process surrogate model is updated by the initial observation records to generate a performance surrogate surface that characterizes the mean and uncertainty of the global performance prediction of the chiller unit under test. By using the performance proxy surface and combining it with the preset thermodynamic constraint boundary of the refrigeration cycle, the extreme value search process is performed on the joint evaluation function consisting of the information collection benefit term and the feasible domain risk penalty term to determine the operating conditions of the next test to be performed. The chiller unit under test is driven to run according to the operating conditions of the next test to be performed and the corresponding performance observation data is collected. The posterior distribution of the performance proxy surface is updated through the collected performance observation data, and the determination and execution of the operating conditions of the next test to be performed are repeated until the total global prediction uncertainty of the performance proxy surface meets the preset termination judgment condition.

[0005] Secondly, embodiments of the present invention provide a working condition self-optimization decision system, including a memory and a processor. The memory stores a computer program that can run on the processor, and the processor executes the program to implement the steps in the above method.

[0006] This invention utilizes a knowledge embedding mapping structure pre-trained from historical test data across different aircraft models to rapidly transfer performance response patterns from historical models to new test models. This allows for the establishment of an initial performance proxy surface with a good global trend under a very small number of initial test conditions, overcoming the problem of blind sampling during the cold start phase of traditional methods. In subsequent testing, a joint evaluation function constructed by fusing thermodynamic constraints of the refrigeration cycle guides the condition search to actively avoid dangerous operating areas and prioritizes exploring regions with high uncertainty and high performance improvement potential, achieving efficient information collection under physical constraints. The performance proxy surface is incrementally updated as testing progresses, and the total global prediction uncertainty is used as the basis for test termination, avoiding invalid testing caused by fixed condition traversal and incomplete characterization due to premature termination. Therefore, this method obtains an accurate description of the full-condition performance of the new model with fewer physical tests, significantly reducing test energy consumption and time, while enhancing the reliability of predictions for unknown operating condition regions and the ability to identify safety risk areas. Attached Figure Description

[0007] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present invention and, together with the specification, serve to explain the technical solutions of the present invention.

[0008] Figure 1 This is a schematic diagram of the software interaction interface provided in an embodiment of the present invention.

[0009] Figure 2This is a schematic diagram illustrating the implementation process of a self-optimization method for chiller performance testing based on cross-model knowledge embedding, provided in an embodiment of the present invention.

[0010] Figure 3 This is a schematic diagram of the hardware entity of a working condition self-optimization decision system provided in an embodiment of the present invention. Detailed Implementation

[0011] This invention provides a self-optimization method for chiller performance testing based on cross-model knowledge embedding. This method can be executed by the processor of a self-optimization decision system. The self-optimization decision system can refer to an industrial control computer or an edge computing terminal. The self-optimization method for chiller performance testing based on cross-model knowledge embedding provided in this invention is deployed in the form of computer program code on the industrial control computer or edge computing terminal, and communicates with the test bench through a graphical user interface to execute the steps of this invention.

[0012] Figure 1The interface shown represents a typical human-computer interaction state during the operation of the method of this invention. The upper left area is the cross-model knowledge embedding mapping structure status display area. This area presents the three-level mapping relationship of the operating condition encoding path, the unit structure encoding path, and the adaptive kernel morphology descriptor fusion path in a flowchart-based manner, and dynamically marks the currently effective structure template channel set and individual specific channel set to reflect the activation status of the cross-model knowledge embedding mapping structure in this test task. The middle right area is the performance proxy surface visualization area. This area displays the global performance prediction mean distribution generated by the Gaussian process proxy model in the form of a heatmap, and uses color depth to represent prediction uncertainty. This surface is dynamically refreshed with the initial observation record and the performance observation data collected in subsequent iterations, corresponding to the steps of generating the performance proxy surface and updating the posterior distribution. The middle left area is the operating condition self-optimization iteration trajectory area. This area marks the measured operating condition points (bubble size represents measured performance value) and recommended trajectory lines in the form of a two-dimensional operating condition space scatter plot, and marks the next operating condition to be tested determined by the current iteration step with a highlighted diamond. This area visually presents the iterative process of extreme value search processing and repeated determination and execution. The bottom left area is the joint evaluation function evolution curve area, which plots the changes of the information collection benefit term, feasible region risk penalty term, and joint evaluation value with each iteration step in real time. This is used to demonstrate the dynamic trade-offs of the joint evaluation function in the extreme value search process and to show the steps of extreme value search processing for the joint evaluation function composed of the information collection benefit term and feasible region risk penalty term. The bottom right area is the system operation log area, which records key events of each working condition recommendation, sampling observation, surrogate model update, and termination condition determination in time series. It particularly highlights the state changes of hidden state tracing and propagation path in the incremental propagation mechanism, which is convenient for operators to verify the execution of the incremental propagation mechanism based on local hidden state tracing. The top status bar displays the current iteration step number, total global prediction uncertainty, current optimal performance value, and termination condition achievement status in real time, which is used to assist in determining whether the termination judgment condition is met. All visual components in the interface are synchronized with the computational kernel of the method. Operators can trigger single-step iterations by clicking on the titles of each module, or the system can automatically execute the repeated steps at preset intervals until the total global prediction uncertainty meets the termination tolerance threshold. At this point, the interface will freeze the performance proxy surface and output the final performance description. This interface serves as both a monitoring panel for the method execution process and an entry point for testers to interact with the self-optimizing algorithm. The following is combined with Figure 2 This invention presents a schematic diagram illustrating the implementation process of a self-optimization method for chiller performance testing based on cross-model knowledge embedding, as provided in an embodiment of the present invention. Figure 2 As shown, the method includes the following steps S100~S500: Step S100: Retrieve the pre-built cross-aircraft type knowledge embedding mapping structure. The cross-aircraft type knowledge embedding mapping structure is obtained by learning and training through meta-tasks from historical multi-aircraft type test sequences. It includes the operating condition encoding path for mapping the operating condition sequence to the operating condition embedding representation, the unit structure encoding path for mapping the unit's inherent attribute information to the structure embedding representation, and the adaptive kernel morphological descriptor fusion path for fusing the operating condition embedding representation and the structure embedding representation and outputting the kernel morphological descriptor set.

[0013] The pre-built cross-model knowledge embedding mapping structure is a general knowledge representation framework extracted from a large amount of historically accumulated performance test data of different models of chillers using machine learning methods. This mapping structure is trained using a meta-task learning strategy, enabling its internal parameters to capture common and transfer patterns across multiple models and multiple test tasks.

[0014] The training process of meta-task learning treats the complete performance test sequence of each historical model as an independent learning task. By iteratively sampling between these tasks, the mapping structure is forced to learn how to quickly build effective performance inference capabilities on unfamiliar models with only a few initial observations. Specifically, this mapping structure consists of three collaborative processing paths. The operating condition encoding path maps the input operating condition sequence to an operating condition embedding representation. Here, the operating condition sequence refers to the combination sequence of external control parameters and internal state parameters of the chiller unit under different operating states, such as the chilled water outlet temperature setpoint sequence, the cooling water inlet temperature sequence, and the compressor load rate sequence. The implementation of the working condition encoding path can be based on a temporal convolutional attention network architecture. This network first extracts multi-scale temporal features from the input working condition sequence through multiple dilated convolutional modules. The dilation factor of the dilated convolutional module increases exponentially with the number of layers, so that the receptive field of the top layer can cover the entire test sequence length. Then, a multi-head self-attention mechanism is used to perform global dependency modeling on the temporal feature map output by the dilated convolution to capture the long-range interaction relationship between different working condition steps. Finally, through a global average pooling layer and a linear projection layer, the variable-length temporal features are compressed into a fixed-dimensional working condition embedding representation.

[0015] The unit structure coding path is used to map the input unit inherent attribute information into a structure embedding representation. The unit inherent attribute information includes static parameters describing the design and construction of the chiller unit, such as compressor type and displacement, heat exchange area and number of passes of evaporator and condenser, form and capacity of throttling mechanism, refrigerant charge, etc. This path can be implemented using a hybrid encoding architecture combining graph neural networks and multilayer perceptrons. First, the refrigeration cycle topology of the unit is abstracted as graph data, where the compressor, evaporator, condenser, and throttling device are nodes, and node attributes are the specifications of the corresponding components. Refrigerant piping and connections are the edges of the graph. The node states are iteratively aggregated and updated through several layers of graph isomorphic networks, ensuring that each node's final representation contains its functional role and neighborhood structure information within the overall cycle. Next, the final representations of all nodes are weighted and summed through a differentiable readout layer. The weights are automatically allocated by an attention mechanism based on the node's impact on global performance, resulting in a fixed-length graph-level structure vector. This graph-level structure vector is then concatenated with the unit's global attributes, such as refrigerant type, total cooling capacity rating, and design energy efficiency ratio, and input into a multilayer fully connected network, ultimately outputting a structural embedding representation.

[0016] The adaptive kernel morphology descriptor fusion path is used to perform cross-domain fusion of operating condition embedding representations and structural embedding representations, and output a set of kernel morphology descriptors. This fusion path does not simply concatenate or add the two embedding representations, but rather uses a decoupling and conditional modulation mechanism to identify the invariant components bound to the refrigeration cycle topology and the model-specific components in the structural embedding representation. These components interact with the smooth evolution trend and local non-stationary fluctuations in the operating condition embedding representation, respectively, thereby generating a set of kernel morphology descriptors that can fully characterize the performance distribution of the current test model in a continuous operating condition space. This set will serve as a core component for subsequently constructing a Gaussian process surrogate model.

[0017] Step S200: Obtain the initial performance observations of the chiller unit under the initial test conditions, and retrieve the inherent attribute information of the chiller unit under test. Using the inherent attribute information as input, generate a structural embedding representation through the unit structure coding path. Using the initial test conditions as input, generate a condition embedding representation through the condition coding path. Perform cross-domain fusion of the structural embedding representation and the condition embedding representation through the adaptive kernel morphological descriptor fusion path, output the initial kernel morphological descriptor set, and use the initial performance observations as the initial observation records.

[0018] In one implementation, step S200 may specifically include the following steps S210 to S260: Step S210: Perform channel attribution and allocation processing on the structural embedding representation. Call the pre-set refrigeration cycle configuration template set in the cross-model knowledge embedding mapping structure. Based on the response tendency of each channel of the structural embedding representation to the refrigeration cycle configuration template, classify the channels into the structural template channel set and the individual specific channel set respectively. Combine the structural template channel set into a structural template representation and combine the individual specific channel set into an individual specific representation. The structural template representation carries the invariant components bound to the refrigeration cycle topology, and the individual specific representation carries the model difference components.

[0019] Channel attribution is a classification operation performed on the channels of the structural embedding representation across various dimensions. The structural embedding representation is generated from the unit structure encoding path; it is a multi-dimensional vector, with each dimension representing a channel. Different channels may capture different aspects of the unit structure. The cross-unit knowledge embedding mapping structure pre-constructs a set of refrigeration cycle configuration templates. This template set is learned from historical multi-unit data during the meta-task learning phase and contains several orthogonal template bases. Each template base corresponds to a typical refrigeration cycle configuration pattern, such as a single-stage compression cycle template, a two-stage compressor economizer cycle template, a flooded evaporator configuration template, and a shell-and-tube condenser configuration template. During channel attribution, the similarity between the response vector of the structural embedding representation in each channel and each refrigeration cycle configuration template is calculated. This similarity can be measured by calculating the projection norm of the response vector onto the template tensor. If the projection norm of a channel onto a specific refrigeration cycle configuration template is significantly higher than that onto other templates and exceeds a preset assignment threshold, it indicates that the channel's activation mode closely matches that configuration template, and the information it carries is strongly correlated with the basic topology of the refrigeration cycle. This channel is then assigned to the structural template channel set. Conversely, if the projection norm of a channel onto all refrigeration cycle configuration templates is relatively uniform and generally lower than the assignment threshold, it indicates that its activation mode does not strictly follow any known typical configuration and tends to reflect the unique differences in component selection, capacity matching, and piping layout of specific models. This channel is then assigned to the individual-specific channel set. After determining each channel individually, all channels belonging to the structural template channel set are spliced ​​together in their original order to form a structural template representation. This representation serves as the conditional modulation signal for subsequent steady-state core substrate generation. Similarly, all channels belonging to the individual-specific channel set are spliced ​​together in their original order to form an individual-specific representation, which carries the personalized characteristics that distinguish this model from the general configuration. Through this channel attribution and allocation process, the structural embedding representation is decoupled into two information streams: invariant components and differential components. This lays the information foundation for subsequent modeling of the steady-state trend of performance and model-specific perturbations.

[0020] Step S220: Perform time-series component separation processing on the working condition embedded representation, extract the smooth-changing main trend representation of the working condition along the working condition sorting direction implied in the working condition embedded representation, and use the extracted residual as the working condition perturbation representation. The main trend representation of the working condition retains the trend of performance evolving smoothly with the working condition, and the working condition perturbation representation retains local non-stationary fluctuations.

[0021] Temporal component separation aims to decompose the performance evolution information contained in the work condition embedding representation into two components with different physical meanings. Although the work condition embedding representation is a fixed-dimensional vector, it is generated by a temporal convolutional attention network that processes the work condition sequence. The activation patterns of each channel within it implicitly encode the order and dependency of performance changes with work conditions along a specific direction. This specific direction is the work condition ranking direction implied in the work condition embedding representation.

[0022] To extract the main trend of smooth evolution along this direction, this step can employ an extraction method combining local weighted regression and total variational denoising. First, the operating condition embedding representation is treated as a multi-dimensional sequence signal along the operating condition sorting direction. Local weighted regression is applied to fit the signal in each channel, using a three-cubic kernel. The bandwidth parameter is adaptively selected based on the span of the entire operating condition space, ensuring the fitted curve maintains the overall trend while being insensitive to local high-frequency fluctuations. The smooth curve obtained after fitting constitutes the main trend representation of the operating conditions. Subsequently, the corresponding channel value of the main trend representation is subtracted from the original operating condition embedding representation channel by channel. The resulting difference sequence is the residual, which is used as the operating condition perturbation representation. The main trend representation captures the physical laws governing the continuous and gradual change of performance with operating conditions, such as the near-linear increase in cooling capacity as the chilled water outlet temperature rises. The operating condition perturbation representation retains local non-stationary fluctuations that cannot be explained by the smooth curve, such as transient performance fluctuations caused by compressor surge or expansion valve oscillation, as well as minor distortions caused by measurement noise. This separation enables subsequent steps to differentiate between the deterministic trend and random perturbation components of the performance.

[0023] Step S230: Using the structural template representation as the conditional modulation signal, retrieve the modulation descriptor set associated with the structural template representation in the cross-model knowledge embedding mapping structure, and perform channel-by-channel readjustment on the main trend representation of the operating condition through the modulation descriptor set, so that the evolution form of the main trend representation of the operating condition is compatible with the cooling cycle configuration, and obtain the steady-state core basis.

[0024] The core idea of ​​the conditional modulation mechanism is to use a structural template representation extracted from the unit's invariant components to adjust the internal feature distribution of the main trend representation of the operating condition, ensuring that the reconstructed trend strictly conforms to the thermodynamic laws that this type of refrigeration cycle should follow. During training, the cross-model knowledge embedding mapping structure learns a set of corresponding modulation descriptors for each typical refrigeration cycle configuration template. These modulation descriptors are stored in the mapping structure as a set of matrices, where each modulation descriptor is a combination of a linear transformation matrix and a bias vector, used to describe a specific modulation operation. When performing channel-by-channel readjustment, the structural template representation is first input into a routing network, which is a lightweight multilayer perceptron. Its output is a set of attention weights, which are applied to all available modulation descriptor sets. Through weighted summation, a customized modulation descriptor sequence best suited to the current unit configuration is dynamically combined from a large pool of descriptors. The sequence length is consistent with the number of channels in the main trend representation of the operating condition. Subsequently, for each channel of the main trend representation of the operating condition, a readjustment operation is performed using the corresponding customized modulation descriptor.

[0025] For the readjustment operation, we can set the scalar value of the i-th channel in the main trend representation of the operating condition as xi, and the corresponding customized modulation descriptor includes a scaling factor αi and a translation factor βi. Then, the readjusted value of this channel is yi = αi * xi + βi. The scaling factor αi and the translation factor βi are both provided by the modulation descriptor combined by the routing network. After channel-by-channel scaling and translation, the amplitude and offset of each channel in the main trend representation of the operating condition are recalibrated, making its evolution highly compatible with the refrigeration cycle configuration of the unit under test. For example, for a configuration using an economizer cycle, the modulation descriptor will enhance the curvature change of the trend curve in the intermediate pressure operating condition region, making it more consistent with the two-stage performance characteristics of the economizer cycle. The readjusted representation is the steady-state core basis, which, as the deterministic skeleton of the performance proxy surface, carries the expected performance form of the unit under ideal, undisturbed conditions.

[0026] Step S240: Perform cross-interaction processing on individual-specific representations and operating condition disturbance representations. Through the preset interaction relationship template in the cross-model knowledge embedding mapping structure, perform signal exchange and component fusion at corresponding positions on individual-specific representations and operating condition disturbance representations to generate disturbance modulation increments that capture the response mode of operating condition disturbances to model-specific behavior.

[0027] The goal of cross-processing is to uncover the unique response patterns that aircraft-specific structural differences exhibit when faced with non-stationary disturbances in operating conditions. Individual-specific representations originate from channels in the structural embedding representation that do not conform to the general configuration, containing personalized design factors specific to the aircraft model; operating condition disturbance representations originate from fluctuations in the operating condition time-series components that cannot be explained by smoothed trends. During the meta-task learning process, the cross-aircraft knowledge embedding mapping structure has learned a series of interaction relationship templates from historical aircraft data. These templates are organized in the form of a key-value pair memory network, where the keys are used to match specific disturbance patterns, and the values ​​store the corresponding specific structural response rules.

[0028] In one implementation, step S240 may specifically include the following steps S241 to S246: Step S241: Segment the individual-specific representation along the channel direction into perturbation response heads. Before segmentation, determine the number of perturbation response heads based on the historical perturbation response diversity indications recorded in the cross-model knowledge embedding mapping structure. This allows each perturbation response head to independently carry a potential perturbation response pattern in a model-specific behavior, resulting in a set of perturbation response heads.

[0029] Individual-specific representations, as multi-dimensional vectors reflecting the design differences of specific aircraft models, may encode the unique response patterns of that model when faced with operating condition disturbances from different physical sources. Examples include the inertial response mode of a large-displacement screw compressor to sudden changes in cooling water temperature, and the rapid heat conduction response mode of a microchannel heat exchanger to flow fluctuations. To decouple these potential, independent model-specific disturbance response modes, the individual-specific representations need to be structurally segmented. The number of segments is not predetermined but dynamically determined by the historical model disturbance response diversity indicator recorded in the cross-model knowledge embedding mapping structure. This historical model disturbance response diversity indicator is a scalar index statistically obtained after the meta-task learning phase, and its value reflects the total number of categories of disturbance response modes caused by structurally specific components among different models in historical multi-model test data.

[0030] Specifically, cluster analysis can be performed on the hidden states of historical models, for example using density-based clustering algorithms. The number of effective clusters obtained from this clustering is an estimate of the diversity of disturbance responses. After determining the number of disturbance response heads, the total number of channels for individual-specific representations is evenly divided into several continuous channel groups. If the number of channels is not divisible, the last group is fine-tuned for alignment. Each group of continuous channels constitutes a disturbance response head, and all disturbance response heads form a disturbance response head set. Thus, each disturbance response head is given the ability to independently carry a potential model-specific disturbance response pattern during initialization, providing a structured prior inductive bias for subsequent refined interaction with the operating condition disturbance representation.

[0031] Step S242: Map the operating condition disturbance representation to generate a disturbance query vector with the same number of disturbance response headers. Each disturbance query vector carries a query description for pattern matching with the corresponding disturbance response header.

[0032] Although the operating condition perturbation representation has the same dimension as the operating condition embedding representation, it carries residual information after removing the smoothing trend, exhibiting stronger locality and high-frequency characteristics. To enable the perturbation response heads to selectively extract perturbation fragments related to their own response patterns, a dedicated query vector needs to be generated for each perturbation response head to actively address the perturbation information of interest. The mapping process is implemented through a two-layer fully connected transform network. The input of this network is the entire operating condition perturbation representation, and the output is a sequence of perturbation query vectors, the number of which matches the set of perturbation response heads. The first layer of this fully connected transform network includes a linear transform followed by a hyperbolic tangent activation function, aiming to compress the high-dimensional perturbation representation into a low-dimensional latent space, enhancing the abstraction of the query. The second layer is a linear projection layer that directly maps the latent space representation into multiple parallel query vectors, with the dimension of each query vector set to be equal to the channel dimension of a single perturbation response head. Through this nonlinear mapping from global disturbance information to local query representation, each generated disturbance query vector contains information as a query description. This query description implicitly specifies what type of disturbance pattern the query vector expects to match in the operating condition disturbance representation, thereby guiding its corresponding disturbance response head to capture the appropriate context.

[0033] Step S243: For each disturbance response header and its corresponding disturbance query vector, perform a mode compatibility comparison between the response modes of each channel in the disturbance response header and the disturbance query vector position by position, and generate the response sensitivity distribution of each part of the disturbance response header to the operating condition disturbance representation based on the comparison results.

[0034] Pattern compatibility alignment is a crucial operation for quantifying the degree of matching between the response patterns represented internally by the perturbation response header and the query description carried by the perturbation query vector. For a given perturbation response header, its shape is a two-dimensional tensor, where rows correspond to channels included in the header in the individual-specific representation, and columns correspond to the internal feature dimensions of these channels. The corresponding perturbation query vector is also a two-dimensional tensor. The alignment operation is performed at the channel granularity of the perturbation response header. Specifically, firstly, the inner product operation is performed between each channel vector of the perturbation response header and the perturbation query vector. The inner product value, after scaling, serves as the initial response strength of that channel to the current query. Subsequently, the initial response strengths of all channels are input into a flexible maximum transfer function layer, which normalizes along the channel dimensions and outputs a probability distribution. Each element in this distribution represents the sensitivity of that channel to the query description. Next, the normalized channel sensitivity is used as a weight to weight the operating condition disturbance representation. This weighting process generates a vector with the same dimension as the operating condition disturbance representation. The value at each position in this vector reflects the degree of attention the disturbance response head pays to different time segments or disturbance components in the operating condition disturbance representation. This position-by-position generated response sensitivity distribution precisely indicates which local segments in the operating condition disturbance are significant and require key extraction, and which can be ignored, for the model-specific behavior handled by the disturbance response head.

[0035] Step S244: Selectively aggregate the operating condition disturbance representation using each group of response-sensitive distributions, extract disturbance fragments compatible with the response mode of each disturbance response head, and arrange all disturbance fragments into a conditional disturbance array according to the disturbance response head.

[0036] Selective aggregation is a soft attention-driven information extraction operation. For each perturbation response head, step S243 has generated a response-sensitive distribution, which has the same sequence length dimension as the condition perturbation representation. The response-sensitive distribution is used as a weight mask and directly multiplied element-wise with the condition perturbation representation. This multiplication operation preserves and enhances the information at positions with high weights under the response-sensitive distribution, while compressing the values ​​at positions with low weights to near zero. After masking and weighting, a weighted summation is performed along the sequence dimension on the weighted condition perturbation representation. The aggregation weights are provided by the response-sensitive distribution itself. The resulting vector is a fixed-dimensional vector, which is the perturbation fragment extracted from the condition perturbation representation. It condenses the perturbation information that the current perturbation response head is interested in and that is highly compatible with its response pattern. After performing this selective aggregation operation on all perturbation response heads, all the resulting perturbation fragment vectors are stacked as row vectors in the original order of the perturbation response heads, forming a two-dimensional conditional perturbation array. The array has the same number of rows as the number of disturbance response headers and the same number of columns as the dimension of each disturbance segment. Each row captures the portion of the disturbance representation related to a specific model-specific behavior.

[0037] Step S245: The conditional perturbation array and the individual-specific representation are fused position by position. During the fusion process, the positions in the individual-specific representation that are not activated by the conditional perturbation array are subject to biased suppression, and the positions in the individual-specific representation that are activated by the conditional perturbation array are enhanced and preserved to obtain the fused perturbation representation.

[0038] In the position-by-position fusion process, the conditional perturbation array is first mapped to a shape identical to the individual-specific representation through a linear transformation layer, resulting in a mapped perturbation array. Subsequently, a gating signal jointly determined by the mapped perturbation array and the individual-specific representation is constructed. This gating signal can be calculated using a single hidden-layer fully connected network, which simultaneously receives the values ​​of the corresponding positions in both the mapped perturbation array and the individual-specific representation, outputting an activation intensity coefficient in the range of zero to one. For each position in the individual-specific representation, if the conditional perturbation array exhibits a high response amplitude at the corresponding position or its vicinity, indicating a strong correlation between that position and the currently captured operational perturbation, the activation intensity coefficient output by the gating network tends to 1, and the value at that position is enhanced and preserved. Conversely, if the conditional perturbation array responds weakly at the corresponding position, indicating that the position does not carry effective information under the current perturbation interaction mode, the activation intensity coefficient output by the gating network tends to 0, and the value at that position is biased and significantly attenuated. The individual-specific representation is multiplied positionally by the corresponding activation intensity coefficient, and then added to the product of the mapped perturbation array and a learnable scaling factor to obtain the final fused perturbation representation. This fused perturbation representation selectively enhances the model-specific components related to the current operating condition perturbation mode, while suppressing irrelevant basis differences, making the representation of the perturbation response mode purer and more prominent.

[0039] Step S246: Apply a time-consistent projection transformation to the fused perturbation representation to map the fused perturbation representation to a time-series arrangement space shared with the steady-state kernel basis, and output the perturbation modulation increment.

[0040] Since the fused perturbation representation originates from the cross-interaction between individual-specific representations and operational condition perturbations, its internal channel arrangement order and numerical distribution range may differ from the steady-state kernel basis. Without alignment, effective fusion in subsequent steps is impossible. The temporal consistency projection transformation aims to establish an isomorphic mapping between the two representation spaces. This projection transformation is implemented through a fully connected layer whose weight matrix is ​​learned during training by a cross-machine knowledge embedding mapping structure, enabling it to learn the affine transformation from the individual-specific space to the steady-state kernel space. The fused perturbation representation, as input, undergoes matrix multiplication and addition operations through this fully connected layer. The output dimension is then adjusted to be completely consistent with the steady-state kernel basis, and the channel arrangement semantics are aligned with the temporal arrangement space of the steady-state kernel basis. Furthermore, to ensure that the magnitude after mapping is consistent with the steady-state kernel basis, a layer normalization operation is cascaded after the fully connected layer. This operation calculates the mean and variance along the channel dimension and standardizes the mapping result, ensuring that the distribution of the output values ​​falls within a controlled range. The output after this projection transformation is the disturbance modulation increment, which accurately measures the modulation offset applied by both model-specific and operating condition disturbances on the expected performance trajectory specified by the steady-state core substrate.

[0041] Step S250: Concatenate the steady-state kernel substrate and the perturbation modulation increment along the channel dimension, and input the concatenation result into the morphological compression mapping. After channel filtering and merging operations, an initial kernel morphological description set with dimensions consistent with the preset kernel morphological description set is obtained. The transformation rules of the morphological compression mapping are determined by the kernel morphological constraints stored in the cross-model knowledge embedding mapping structure.

[0042] The steady-state kernel basis represents the deterministic framework for smooth performance evolution, while the perturbation modulation increment represents the aircraft-specific deviation response to perturbations. Both characterize the performance profile of the tested aircraft from different perspectives. Cascading them along the channel dimension means directly concatenating their channels while maintaining consistency in other dimensions, forming a joint representation tensor with the sum of their channels. This joint representation tensor is then subjected to input morphological compression mapping. Morphological compression mapping is a lightweight compression network based on alternating depthwise separable convolutions and pointwise convolutions. Its design goal is to merge and filter redundant or highly correlated information channels in the joint representation tensor, compressing them into a compact representation with dimensions strictly consistent with the preset set of kernel morphological descriptors. This compression network integrates kernel morphological constraints stored in the cross-aircraft knowledge embedding mapping structure. These constraints manifest as spectral norm regularization terms for the convolution kernel weights and penalty terms for correlations between output channels in each compression layer. They ensure that the compression process does not destroy the mathematical properties of the kernel descriptors, such as symmetry, positive definiteness, and smoothness. The first stage of the compressed network is a pointwise convolutional layer, which initially compresses the number of channels to an intermediate dimension through linear combination across channels. This is followed by several residual blocks composed of depthwise separable convolutions, batch normalization, and corrected linear activation units, further reorganizing features while preserving spatial information. Finally, a pointwise convolutional layer without activation precisely maps the number of channels to a preset dimension, and a global response linear normalization layer calibrates the range of output values. The compressed output is the initial set of kernel morphology descriptors, a compact multidimensional tensor that can be directly used as the kernel function parameterization basis for the Gaussian process surrogate model in subsequent steps.

[0043] In one implementation, step S250 may specifically include the following steps S251 to S256: Step S251: Apply a numerical range equalization operation to the steady-state kernel basis, adjust the numerical distribution of each channel to make it fall into a consistent distribution range, and obtain an equalized steady-state basis. The distribution range is determined by the preset steady-state basis distribution specification in the cross-model knowledge embedding mapping structure.

[0044] During the generation of the steady-state kernel basis, due to the inherent differences in physical dimensions and sensitivity of the performance trend information carried by its channels, the mean and variance of the numerical distribution of different channels may vary significantly. If directly fused with subsequent representations, channels with higher numerical magnitudes will dominate in the compression mapping, weakening the balance of information representation. The numerical range balancing operation addresses this problem. It does not change the relative structure within each channel of the steady-state kernel basis, but only performs an independent affine transformation on each channel to adjust its numerical distribution to a consistent range. The cross-model knowledge embedding mapping structure pre-defines a set of steady-state basis distribution specifications. This specification is derived at the end of the meta-task learning by statistically analyzing the median and half-quartile of the distribution of each channel of the steady-state kernel basis generated under ideal conditions for all historical models, and is stored as a set of target mean vectors and target standard deviation vectors. During the balancing operation, for each channel of the steady-state kernel basis, its current mean and standard deviation are first calculated, and then the transformation is applied: the new value equals (old value minus current mean) divided by the current standard deviation, then multiplied by the target standard deviation, and finally added to the target mean. After this channel-by-channel standardization and rescaling operation, the numerical distribution of all channels is pulled into a consistent target distribution range, ensuring that each channel is treated equally during subsequent cascading and compression, resulting in a balanced steady-state basis.

[0045] Step S252: Perform boundary constraint processing on the perturbation modulation increment. Based on cross-model knowledge, embed the perturbation tolerance boundary stored in the mapping structure to limit the response amplitude of the perturbation modulation increment that exceeds the boundary to the boundary edge, and obtain the corrected perturbation increment.

[0046] The perturbation modulation increment captures performance deviations caused by local non-stationary fluctuations and model-specific characteristics. Because the generation process involves nonlinear interactions, its output response may exhibit abnormally large peaks in rare cases. If these peaks, exceeding the normal physical tolerance range, are directly fed into subsequent processing, they may cause unreasonable distortions in the proxy surface. The cross-model knowledge embedding mapping structure stores perturbation tolerance limits. These limits are a set of absolute numerical thresholds corresponding to each channel of the perturbation modulation increment, statistically determined based on the maximum perturbation modulation amplitude actually observed in numerous tests of historical models. Limit constraint processing employs a truncation function. For each element in the perturbation modulation increment tensor, if its value is greater than the upper limit of the perturbation tolerance for the corresponding channel, its value is hard-set to the upper tolerance limit; if its value is less than the lower limit of the perturbation tolerance for the corresponding channel, its value is set to the lower tolerance limit; if its value is within the limits, its original value remains unchanged. While this hard truncation method introduces slight discontinuities at the boundaries, it effectively prevents individual outliers from causing deviations in the overall performance inference process. The result of this processing is the corrected disturbance increment, which retains the effective information of the disturbance response while limiting its extreme amplitude.

[0047] Step S253: Connect the equilibrium steady-state basis and the correction perturbation increment along the channel direction, and insert a transition bridging channel at the connection point. The value of the transition bridging channel is jointly determined by the adjacent channels of the equilibrium steady-state basis and the correction perturbation increment, which is used to mitigate the information gap between the two parts and obtain the joint representation tensor.

[0048] Directly concatenating the equilibrium steady-state basis and the correction perturbation increment along the channels creates a sharp information boundary at their intersection. The smooth, deterministic characteristics of the equilibrium steady-state basis abruptly change the volatility of the correction perturbation increment, hindering smoothing of subsequent compression mapping. Therefore, a small number of transition bridging channels are inserted at the connection. The number of transition bridging channels can be fixed at a small value, such as 5% of the total. The value of each transition bridging channel is generated by a fusion formula inspired by bilinear interpolation: for each transition bridging channel i, its value is (1-λi) multiplied by the last channel of the equilibrium steady-state basis, plus λi multiplied by the first channel of the correction perturbation increment. Here, λi is an interpolation coefficient linearly increasing from 0 to 1, assigned to different transition bridging channels. In this way, the initial transition channel content is close to the equilibrium steady-state basis, the final transition channel content is close to the correction perturbation increment, and the intermediate content is smoothly transitioned. These inserted channels have no physical counterpart, but their smoothing effect ensures information continuity in the entire cascaded tensor along the channel direction, mitigating the information gap between the two main parts. After insertion, the equilibrium steady-state basis, the transition bridging channel sequence, and the correction perturbation increment are connected end to end along the channel dimension to form a complete joint representation tensor.

[0049] Step S254: Compress and map the input form of the joint representation tensor, and reduce the number of channels step by step through multiple compression units. Each compression unit contains staggered one-dimensional channel convolution and pointwise nonlinear mapping, and records the distribution pattern of active channels after each compression.

[0050] The morphological compression mapping is implemented using a deep network consisting of stacked compression units. Each compression unit contains two core operations. The first operation is a one-dimensional channel convolution, for example, using a one-dimensional convolutional layer with a kernel size of 3, sliding along the channel dimension to linearly combine the information from three adjacent channels, reducing redundancy between channels. Because it is a one-dimensional channel convolution, it only performs convolution in the channel direction, without changing the size of other dimensions, effectively reducing the number of parameters to be learned and enhancing the ability to extract local channel features. The second operation is a pointwise nonlinear mapping, consisting of a fully connected layer followed by a hyperbolic tangent activation function. The fully connected layer performs the same linear transformation and nonlinear activation on the cross-channel vector at each independent position, achieving pointwise feature transformation and enhancing the nonlinear expressive power of the representation. One-dimensional channel convolution and pointwise nonlinear mapping are alternately arranged within this compression unit, i.e., channel convolution first, then pointwise mapping, then channel convolution again, then pointwise mapping again. After processing by this compression unit, the number of channels is compressed to about two-thirds of the original. After each level of compression, the output representation of that level undergoes active channel filtering. Specifically, the variance of the activation value for each channel is calculated, and the channels with the highest variance are marked as active channels for that level. The distribution pattern index of these active channels in the original channel arrangement is recorded. This index records which channels are determined to be active channels with strong information carrying capacity. This step-by-step recording provides an indexing mechanism for achieving the interpretability and traceability of the final kernel morphological description.

[0051] Step S255: At the output of each compression unit, apply symmetry adjustment to the current compression representation, adjust the interrelation structure between channels, so that the correlation morphology matrix that can be constructed based on the compression representation falls into the symmetric positive definite domain, and obtain the symmetric compression representation.

[0052] The kernel morphology descriptor set is used to define the correlation measurement rules between different operating points in a Gaussian process. These correlation measurement rules require that the corresponding correlation morphology matrix be a symmetric and positive definite matrix; otherwise, numerical computation instability will occur. Therefore, after each compression unit output, the current compressed representation is symmetricized. This symmetry-based projection algorithm can be used. Let the current compressed representation be a two-dimensional matrix, where rows can be considered as representations of different virtual points and columns as channels. First, the Gram matrix of this matrix is ​​calculated, which is the product of the matrix and its transpose. This Gram matrix represents the cross-correlation structure between channels. Theoretically, the final correlation morphology matrix should be derived from this Gram matrix, so it needs to be forced into a symmetric positive definite region. By performing eigenvalue decomposition, the Gram matrix is ​​decomposed into an eigenvector matrix and an eigenvalue diagonal matrix. Then, all negative eigenvalues ​​and eigenvalues ​​less than a small positive threshold are forcibly replaced with the threshold itself, while positive eigenvalues ​​greater than the threshold remain unchanged. Afterward, a new, positive definite Gram matrix is ​​reconstructed using the original eigenvector matrix and the corrected eigenvalue diagonal matrix. Finally, we find a matrix such that the product of this matrix and its transpose approximates the modified Gram matrix as closely as possible. This search can be achieved through a convergent iterative process, with each step performing a least-squares approximation under orthogonal constraints. The matrix obtained after convergence is the symmetric compressed representation after symmetry simplification, ensuring that any correlational form matrix constructed based on it naturally lies within the symmetric positive definite region.

[0053] Step S256: Output the final first-level symmetric compression representation as the initial kernel morphological descriptor set, and at the same time save the distribution pattern index of the active channels in the last layer of the morphological compression mapping.

[0054] After multi-level compression, the final stage outputs a symmetric compressed representation whose channel dimensions perfectly match the dimensions of the preset kernel morphological descriptor set. This symmetric compressed representation, refined through multi-level information extraction and channel filtering, highly condenses the performance morphological knowledge jointly expressed by the steady-state kernel basis and perturbation modulation increments. It also naturally possesses the mathematical properties of a Gaussian process kernel function parameter, and is therefore directly output as the initial kernel morphological descriptor set. Simultaneously, the morphological compression mapping records the distribution pattern index of active channels at each level. In the final layer, this index indicates which dimensions have higher information activity and contribute more to describing the performance morphology within the final kernel morphological descriptor set. Saving this final layer's active channel distribution pattern index along with the initial kernel morphological descriptor set provides fine-grained channel-level guidance for subsequent model diagnosis, kernel morphological interpretation, and knowledge solidification operations after testing.

[0055] Step S260: Write the initial performance observations, the condition embedding representations corresponding to the initial test conditions, and the initial kernel morphology description set into the observation cache to obtain the initial observation record containing the mapping relationship between the condition location identifier and the performance observations.

[0056] The observation cache is a memory-based data storage structure used to accumulate observation information for all executed operating points throughout the entire testing process. For the initial test involved in this step, the following are written to the observation cache as a triple: the initial performance observations collected during the test (i.e., the specific values ​​of cooling capacity, compressor power consumption, and coefficient of performance); the operating condition embedding representation generated after inputting the initial test operating condition into the operating condition encoding path (this embedding representation serves as the unique location identifier of the operating point in the embedding space); and the initial kernel morphology description set generated in the preceding steps. Each record in the cache is indexed by the operating condition location identifier, which is the operating condition embedding representation, while the performance observation mapping relationship refers to the performance observation value corresponding to that operating point. After writing, the current state of the cache constitutes the initial observation record. This record, representing all the knowledge of the tested unit's performance to date, will be used to initialize the Gaussian process surrogate model, providing prior conditions for its subsequent updates.

[0057] Step S300: Construct a Gaussian process surrogate model based on the initial set of kernel morphological descriptors, update the posterior distribution of the Gaussian process surrogate model through the initial observation records, and generate a performance surrogate surface that characterizes the mean and uncertainty of the global performance prediction of the chiller unit under test.

[0058] The Gaussian process surrogate model is a nonparametric Bayesian regression model. Its core lies in using a kernel function to define the correlation between any two points in the operating space and inferring the function value distribution of the unobserved region based on observed data points. An initial set of kernel morphological descriptors is used as the parameter source for defining the kernel function of this Gaussian process surrogate model. Specifically, let x and x' be the coordinate vectors of two operating points. The value of their kernel function k(x,x') is calculated by a metric function determined by the initial set of kernel morphological descriptors. This metric function can be considered as an embedding kernel. First, x and x' are each mapped to the latent space through an up-dimensional embedding network parameterized by the initial set of kernel morphological descriptors. Then, the inner product of the two vectors in the latent space is calculated, followed by a nonlinear exponential squared kernel transformation to obtain the final kernel function value. The initial observation records contain operating point identifiers and their corresponding performance observations. Using these initial observation records as training data, and based on Bayes' theorem, the prior Gaussian process is combined with the likelihood of the observed data to calculate the posterior distribution. The posterior distribution update calculation can be performed in function space using a stable square root unscented Kalman filter or its variants. This ultimately calculates the performance prediction mean for each node on the global operating condition grid, representing the expected estimate of the node's performance coefficients, and simultaneously calculates the prediction uncertainty, typically expressed as the prediction variance, which is a measure of the reliability of the estimate. These two quantities together depict a surface on the global operating condition grid, known as the performance proxy surface. This surface not only indicates the most likely performance of the tested unit under different operating conditions but also clearly identifies areas where performance inferences lack observational data support and have significant uncertainty.

[0059] In one implementation, step S300 may specifically include the following steps S310 to S360: Step S310: Configure a set of induced positions in the operating space of the chiller unit to be tested, and assign a hidden state description to each induced position. The initial value of the hidden state description is determined by the prior morphology defined by the initial kernel morphology description set, and the initial association between each induced position is determined by the association measurement rule given by the initial kernel morphology description set.

[0060] To enable the Gaussian process surrogate model to efficiently handle the rapidly increasing computational complexity with the number of observations, a sparse Gaussian process framework based on induced points is adopted. Configuring induced positions refers to placing a set of virtual reference points—far fewer than the total number of nodes in the global grid—within the continuous operating space of the entire chiller under test, either fixedly or adaptively, according to a certain strategy. These virtual reference points are the induced positions. Initial point placement can employ Latin hypercube sampling or Haldane quasi-random sequence sampling strategies to cover the entire operating space distribution as comprehensively as possible with the fewest possible points. A hidden state description is assigned to each induced position. This hidden state description is a low-dimensional vector containing sufficient statistics for the mean vector and covariance matrix; the initial value setting of this description is crucial. The initial values ​​are determined by the prior morphology defined by the initial kernel morphological descriptor set. Specifically, the working condition coordinates of the induced positions are input into the embedding network parameterized by the initial kernel morphological descriptor set to obtain a prior hidden state representation. This representation is then converted into sufficient statistical initial values ​​for the hidden state description through a linear readout layer. The mean vector is the first part of the output of the linear readout layer, and the initial form of the covariance matrix is ​​a diagonal matrix. The diagonal elements are obtained by exponentializing the second part of the output of the linear readout layer to ensure positive definiteness. The initial correlation between each induced position is determined by the correlation metric rule given by the initial kernel morphological descriptor set. Substituting the working condition coordinates of any two induced positions into the aforementioned embedding kernel calculation, the resulting kernel function value is filled into an initial prior correlation matrix, serving as the basis for defining the conditionally independent structure among these hidden state descriptions.

[0061] Step S320: Using the working condition location identifier and performance observation in the initial observation record as conditions, generate the response likelihood relationship between the hidden layer state description of each induced location and the observation record, and combine the response likelihood relationship with the prior association morphology derived from the initial set of kernel morphology descriptions to obtain the structured joint distribution.

[0062] The response likelihood relation defines the probability distribution of observed performance values ​​at each working point location in the initial observation record for any induced location, given its hidden state description. This likelihood is typically assumed to be Gaussian, with its mean output by a decoding network. The decoding network takes the hidden state description of the induced location and the target observed working point location identifier as inputs, calculates the interaction between the hidden state and the working point embedding, and outputs the predicted mean performance. The variance of the likelihood is a quantity controlled by an independent noise parameter, initially set to a small positive constant. For each induced location, all observation data from the initial observation record are substituted into this Gaussian likelihood to obtain a series of likelihood functions. The prior correlation morphology is the covariance matrix composed of the kernel function values ​​between each pair of induced locations calculated using the embedding kernel from the initial kernel morphology description set. It specifies the multivariate Gaussian prior distribution of the hidden state description at each induced location when no observation data has been observed. Multiplying the response likelihood relation with the prior association pattern, i.e., following the chain rule of probabilistic graphical models, constructs a complete joint probability distribution containing all induced location hidden state descriptions and all observation performance values, i.e., a structured joint distribution. The structured nature of this distribution is reflected in the fact that the association between the induced location hidden state descriptions strictly follows a conditionally independent graph defined by the kernel function, while the observations are conditionally independent of the hidden state, allowing subsequent inferences to fully utilize sparsity.

[0063] Step S330: Under the premise of maintaining the independent structure of the hidden state descriptions at each induced position as specified by the prior association morphology, perform variational inference on the structured joint distribution, iteratively adjust the sufficient statistics of each hidden state description, so that the lower bound of variational evidence increases in round by round, and obtain the variational posterior hidden state set.

[0064] In one implementation, step S330 may specifically include the following steps S331 to S335: Step S331: Initialize the variational distribution configuration of the hidden layer state description at each induced location. Take the observation samples that are closest to each induced location in terms of working space distance from the initial observation record as local evidence. Use the local evidence to locally warm up the variational distribution configuration so that the variational distribution is initially biased towards local observations.

[0065] The initialization quality of variational inference directly affects the convergence speed and the quality of the final solution. Therefore, instead of initializing the variational distribution of all induced locations to a uniform prior distribution, a warm-up strategy based on local nearest neighbors is implemented. First, for each induced location, the Euclidean distance between its working condition coordinates and the working condition location identifiers of all observation samples in the initial observation record is calculated, and the K nearest observation samples are selected, where K can be between 3 and 5. The performance observations of these K local nearest neighbor observation samples are used as local evidence for that induced location. The initialization process of the variational distribution configuration is as follows: the variational mean vector is initialized as the weighted average of these K performance observations, with the weights proportional to the reciprocal of the Euclidean distance from each observation sample to the induced location; the variational covariance matrix is ​​initialized as a diagonal matrix, with the diagonal elements first set to the sample variance of these K performance observations, and then multiplied by a scaling factor slightly greater than 1 to reflect additional uncertainty. Through this local preheating, the variational distribution at each induced location carries the observation information of its neighboring region before the global inference loop is started, which makes the variational distribution initially biased towards local observations and greatly accelerates the convergence of subsequent global consistent inference.

[0066] Step S332: Calculate the mutual information between the hidden state description of each induced position and the global observation sample set under the current variational distribution. Sort the influence of each induced position according to the mutual information and select the induced positions with the highest mutual information to form the active induced set of the current iteration round.

[0067] Not all induced locations contribute equally to the current inference task. Some induced locations located in sparsely observed regions have hidden state descriptions primarily determined by priors, with minimal updates and weak influence on global inference. To identify the truly crucial induced locations in each iteration, the mutual information between the hidden state descriptions of each induced location and the global observation sample set is calculated. Under the current variational distribution, this mutual information can be approximated by calculating the difference between the entropy of the variational distribution and the variational conditional entropy of the observed data given the hidden state description. Higher mutual information indicates that the hidden state description of the induced location is more critical to interpreting the observed data. After calculating the mutual information of all induced locations, they are sorted in descending order, and a selection ratio is set, for example, selecting the top 50% of induced locations with high mutual information, marking them as the active induced set for the current iteration. This influence-based selection mechanism based on mutual information allows the computational resources for variational inference to focus on the most information-rich regions that require key updates.

[0068] Step S333: For the induction positions in the active induction set, calculate the gradient signal of the sufficient statistic describing the hidden state with respect to the lower bound of the variational evidence. For the induction positions in the inactive induction set, keep the sufficient statistic unchanged to maintain the stability of the global structure.

[0069] In step S333, the gradient signal is the key guide that drives the sufficient statistics of the variational distribution to update in the direction of maximizing the lower bound of variational evidence. In this iteration, the computational operation is limited to the active inducement set. For each inducement position in the active inducement set, using automatic differentiation, the gradient of the lower bound of variational evidence with respect to the elements of its variational mean vector and variational covariance matrix is ​​calculated at the current sufficient statistics of variation. This gradient signal indicates in which direction these sufficient statistics should be fine-tuned to raise the lower bound. For inducement positions assigned to the inactive inducement set, the sufficient statistics of their variational distribution are frozen in this iteration, maintaining their original values. This selective update strategy significantly reduces the computational cost of a single iteration. Furthermore, by preserving the state of most inactive positions, it can serve as a stable anchor point, maintaining the relative stability of the global structure under the constraints of prior correlation morphology, preventing the overall covariance structure from being torn apart by drastic local updates, and thus playing a role in smoothing inference.

[0070] Step S334: Perturb and update the sufficient statistics of the active induced positions using the gradient signal, and apply a structure preservation constraint after the update to ensure that the relative distance between the updated hidden state descriptions of the active induced positions is consistent with the similarity specified by the prior association morphology.

[0071] In step S334, the gradient signal calculated in step S333 is used to perform a perturbation update on the sufficient statistics of the active induced positions. The update employs an adaptive step-size stochastic gradient ascent method or its momentum variant, such as the Adam optimizer, using the exponential moving average of the gradient as the first and second moment estimates. After correcting the gradient bias, the parameter update is calculated and applied to the sufficient statistics. To maintain the structural integrity of the Gaussian process surrogate model during the update process, a structure-preserving constraint needs to be applied to the relative relationships between the hidden state descriptions of the active induced positions after the update. This constraint is achieved by adding a regularization term to the loss function, i.e., the lower bound of the negative variational evidence. This regularization term penalizes the deviation between the cosine similarity between any two active induced positions in the current variational mean vector and the normalized kernel function value of that position in the prior association pattern. By encoding this structure into the objective function, the gradient signal itself contains directional guidance to maintain the structure. While the updated active induced position hidden state descriptions tend towards the data, their relative distances and similarity patterns remain consistent with the prior similarity patterns characterized by the initial set of kernel morphology descriptions, ensuring the consistency of the inference results in both physical and statistical senses.

[0072] Step S335: Check the improvement of the variational evidence lower bound after the update. When the improvement is no longer significant, prune the inactive induced positions from the subsequent calculations and fix the set of induced positions retained after pruning as the new induced layout. Output the final variational posterior hidden state set, where the sufficient statistics describing each hidden state include the inverse representation of the local predicted mean vector and the local correlation matrix at that induced position.

[0073] Variational inference is a progressively converging process. After each iteration, the relative improvement of the variational evidence lower bound between the current and previous iterations is calculated. If this relative improvement remains below a preset convergence tolerance threshold for multiple iterations (e.g., three), the variational inference is considered converged. At this point, induction positions that have long been inactive induction sets and contribute negligibly to the inference are removed from the entire induction point layout—a pruning operation. The pruned set of induction positions retains only the points with the highest information content describing the performance proxy surface, and is fixed into a new, more streamlined induction layout. In this layout, each retained induction position possesses a converged variational posterior hidden state description, and the variational mean vector in its sufficient statistics is the local prediction mean vector; the inverse of its variational covariance matrix, i.e., the inverse representation of the local correlation matrix, is also calculated and stored. The set of these sufficient statistics constitutes the final variational posterior hidden state set, which sparsely but highly representatively summarizes the information passed from all observation data to the Gaussian process surrogate model.

[0074] Step S340: Using the variational posterior hidden state set as the source, and the correlation metric rule implied in the prior correlation form as the diffusion conductance tensor, the manifold is diffused into the continuous working space, propagating the hidden state at the discrete induced position to every node of the global working grid, generating a global hidden state manifold.

[0075] The variational posterior hidden state set only provides state information at a finite number of discrete points. To generate a global performance proxy surface, this discrete information needs to be made continuous. This propagation process is implemented using a manifold diffusion method based on kernel regression. Discrete induced positions are considered as source points, each carrying its variational posterior hidden state description; each node on the global load cell grid is considered as a target point. The inverse representation of the local correlation matrix provided by the variational posterior hidden state set is used as a constraint, and the correlation metric rule implied in the prior correlation morphology is used as the diffusion propagation tensor. The diffusion propagation tensor is defined as the result of evaluating the embedded kernel operator between the target and source points, determined by the initial kernel morphology description set. For each node on the global load cell grid, its hidden state is obtained by weighted average interpolation of the hidden state descriptions of all source points, and its weight is determined by the diffusion propagation tensor: the weight is equal to the normalized result of the kernel function values ​​between the target and each source point. Through this kernel-based smooth interpolation diffusion, the local predicted mean vectors and uncertainty information at all discrete induced locations are propagated to every node of the global load cell grid, forming a vector field that smoothly varies across the entire continuous load space, i.e., the global hidden state manifold. This manifold fully characterizes the surrogate model's understanding of the potential performance state at each load location under posterior inference.

[0076] Step S350: On the global hidden layer state manifold, calculate the performance prediction mean node by node, and calculate the prediction uncertainty based on the local diffusion conductivity tensor trace of the node to obtain the performance surrogate surface indexed by the operating condition coordinates.

[0077] Each node on the global hidden state manifold possesses an interpolated hidden state vector, which is a compressed representation of the performance-related latent states at that node. To obtain the physically meaningful performance prediction mean, the performance coefficients need to be decoded from this hidden state vector. Decoding is accomplished using a linear regression readout, whose parameters are part of a Gaussian process surrogate model and are learned along with variational inference during initialization. The performance prediction mean for each node is obtained by taking the inner product of the hidden state vector of each node with the weights of the readout and adding a bias. Meanwhile, the prediction uncertainty depends not only on the uncertainty of the hidden state estimation but also on the reliability of the node's location in the condition space itself. The prediction uncertainty is calculated based on the local diffusion conductivity tensor trace of the node. Here, the tensor trace specifically refers to the trace of the square matrix formed by the inner product of the interpolated weight vector and its transpose, and its value reflects the sum of the intensity of information diffused from neighboring source points to the node. The prediction uncertainty is ultimately expressed as a sum consisting of the prior noise variance, the variational posterior covariance projection under interpolation, and the reciprocal of the diffusion conductivity tensor trace. By storing the performance prediction mean and prediction uncertainty of all nodes on the global load cell grid in a one-to-one correspondence, a data structure indexed by load cell coordinates is formed. This data structure spans two surfaces in three-dimensional space, representing the performance prediction mean and prediction uncertainty respectively, collectively referred to as the performance surrogate surface.

[0078] Step S360: Traverse the performance proxy surface, mark the node prediction uncertainty corresponding to the region outside the boundary that intersects with the thermodynamic constraint boundary of the refrigeration cycle as invalid, and remove the marked invalid regions from the search domain of the subsequent joint evaluation function.

[0079] The thermodynamic constraint boundary of the refrigeration cycle is a closed surface defined by the first and second laws of thermodynamics, representing the physically achievable operating limits of the chiller unit. Outside this boundary lies a physically inaccessible operating region, where predictions from the performance proxy surface are meaningless extrapolations. This step examines each node on the performance proxy surface. If its coordinates lie outside the thermodynamic constraint boundary, its prediction uncertainty is uniformly marked as a predefined maximum invalid value, such as the maximum value of a floating-point number, regardless of the node's actual prediction uncertainty. This marking operation essentially blacklists these invalid regions. Subsequent steps, when performing extreme value searches to determine the next test condition, first check the search domain. Any region marked as invalid will have its corresponding joint evaluation value set to a minimum or excluded from comparison, thus strictly eliminating it from the search domain. This ensures that the self-optimization process of the test condition always occurs within the physically feasible operating range of the unit, preventing the determination of unexecutable test points due to mathematical inference overflow.

[0080] Step S400: Using the performance proxy surface and combined with the preset thermodynamic constraint boundary of the refrigeration cycle, perform extreme value search processing on the joint evaluation function consisting of the information acquisition benefit term and the feasible region risk penalty term to determine the operating conditions of the next test to be performed. The information acquisition benefit term is determined based on the improvement magnitude of the prediction uncertainty and the prediction mean, and the feasible region risk penalty term is determined based on the degree of violation of the thermodynamic constraint boundary.

[0081] In one implementation, step S400 may specifically include the following steps S410 to S460: Step S410: Discretize the thermodynamic constraint boundary of the refrigeration cycle into a set of boundary points, configure a repulsive response core for each boundary point, and the intensity of each repulsive response core decreases with the distance from the corresponding boundary point. The superposition of all repulsive response cores forms a feasible region risk penalty field on the global operating condition grid. The distribution intensity of the feasible region risk penalty field is determined by the proximity to the constraint boundary.

[0082] The preset thermodynamic constraint boundary of the refrigeration cycle is a continuous surface defined by a mathematical expression, which needs to be transformed into a field that can act on discrete grid points. First, samples are taken on the boundary surface at equal area or equal angle intervals to obtain a large number of dense discrete boundary points, forming a boundary point set. Then, a repulsive response kernel is configured for each boundary point. This kernel function can be either a Laplace kernel function or a Cauchy kernel function centered at that boundary point. Its function value monotonically decreases with the Euclidean distance between the evaluation point and the boundary point, and the decay rate is controlled by a preset bandwidth parameter, which determines the deterrent range of the constraint boundary. For each node on the global operating condition grid, the repulsive response kernel values ​​of all boundary points to it are accumulated and calculated. The sum is the feasible region risk penalty field intensity at that node. Obviously, the closer the grid node is to the constraint boundary, the more high-intensity kernel responses it receives from more boundary points, and the higher its feasible region risk penalty field value; while for nodes deep in the feasible region and far from the boundary, the penalty field value tends to 0.

[0083] Step S420: On the performance proxy surface, extract the predicted mean of each point in the global working condition grid, compare it with the currently discovered global optimal predicted mean, retain the positive difference part, and combine the positive difference part with the prediction uncertainty of the point to generate the original information gain signal.

[0084] Information gain is a preliminary quantification of the profit that can be gained from testing a point. First, the entire grid of test cases is traversed to obtain the average performance prediction for each point. This average can be the average predicted performance coefficient. Simultaneously, a currently known globally optimal average prediction is maintained and continuously updated. This value is the optimal performance coefficient estimate after observation correction among all tested test points to date. For each grid point, the difference between its average performance prediction and the current globally optimal average prediction is calculated. If the difference is greater than 0, it indicates that the point has the potential to discover better performance, and this positive difference is retained. If the difference is less than or equal to 0, it indicates that the predicted performance of the point is not better than the currently known best, and its positive difference is set to 0. Then, this positive difference is combined with the prediction uncertainty of the point, which can be done by multiplication, to generate the original information gain signal. This combination means that the highest original information gain should only be given to a point when it not only has higher predictive performance but also has significant uncertainty, meaning that its authenticity has not been fully verified. Conversely, if the predictive performance is high but the uncertainty is low, it means that the area has been fully explored and confirmed, and its gain should be reduced.

[0085] Step S430: Convolve the original information gain signal with a preset gain pattern modulation function to obtain the information acquisition gain distribution. The gain pattern modulation function has the pattern characteristic of gain attenuation near the current sampled area and gain enhancement in the moderate distance area.

[0086] In one implementation, step S430 may specifically include the following steps S431 to S436: Step S431: Treat the distribution of the original information gain signal on the global working condition grid as a two-dimensional gain distribution map and extract its spatial structure.

[0087] Each node on the global load cell grid possesses an original information gain signal value. These values ​​collectively constitute a scalar field defined on a two-dimensional load cell coordinate network, i.e., a two-dimensional gain distribution map. To understand the essence of its spatial structure and extract its spatial structural features, this can be accomplished by calculating the gradient magnitude field and the Hersey matrix eigenvalue field of this two-dimensional gain distribution map. The gradient magnitude field reveals the degree of drastic local changes in the gain signal, while the Hersey matrix eigenvalue field indicates whether the point is in a local maximum, minimum, or saddle point. Statistical analysis of these geometric features yields macroscopic parameters such as the average gradient and peak density. These parameters reflect the spatial clustering and dispersion of the original information gain signal, providing a basis for the next step of precisely setting the shape parameters of the modulation function.

[0088] Step S432: Construct the profit pattern modulation function as a ring distribution kernel with a central suppression region and a surrounding enhancement band. The radius of the central suppression region of the ring distribution kernel is determined by the minimum envelope radius of the currently sampled operating point in the operating space, and the radius of the surrounding enhancement band is determined by the average amount of the local correlation length in the performance proxy surface.

[0089] The profit pattern modulation function is constructed as a binary ring-distributed kernel function with a well-defined geometric partition in space.

[0090] When calculating the radius of the central suppression region, the coordinates of all currently sampled operating points in the operating space can be obtained first. The minimum envelope hypersphere radius of this coordinate point set is then calculated, which is the Euclidean distance from the geometric center of the point set to the furthest sampled point. This minimum envelope hypersphere radius is multiplied by a suppression expansion coefficient to obtain the radius of the central suppression region. This suppression expansion coefficient is provided by a pre-defined exploration tendency parameter table in the cross-model knowledge embedding mapping structure. This exploration tendency parameter table is obtained through joint optimization of grid search and cross-validation during the meta-task learning phase. Its specific value is determined by a lookup table indexed by the total global prediction uncertainty of the current performance proxy surface and the average number of termination iterations in historical model tests. When the total global prediction uncertainty is high and the historical average number of termination steps is large, a smaller suppression expansion coefficient is used to encourage more aggressive exploration; conversely, a larger value is used for conservative contraction.

[0091] When determining the radius of the surrounding enhancement band, we can first extract the implicit length scale parameter from the set of current kernel morphological descriptors associated with the performance proxy surface. Then, we perform spectral decomposition on each channel of the kernel morphological descriptor set, taking the square root of the principal eigenvalue of each channel as the estimate of the local correlation length. Next, we calculate the arithmetic mean of the local correlation length estimates for all channels to obtain the average local correlation length. Multiplying this average by an enhancement scaling factor gives us the outer radius of the surrounding enhancement band. This enhancement scaling factor is also provided by the exploration tendency parameter table, and its indexing method is consistent with the suppression expansion coefficient. Furthermore, it exhibits an inverse relationship with the suppression expansion coefficient; that is, when the suppression expansion coefficient takes a larger value, the enhancement scaling factor takes a smaller value to maintain a dynamic balance between exploration and conservatism in the geometric proportion structure between the central suppression region and the surrounding enhancement band. The resulting ring-shaped distribution kernel fully characterizes the desired gain modulation spatial morphology.

[0092] Step S433: Perform a spatial convolution operation between the two-dimensional revenue distribution map and the annular distribution kernel. This reduces the original information revenue in densely sampled working point regions and amplifies it in sparsely sampled working point regions that maintain a suitable distance from the sampled areas, thus obtaining an enhanced revenue distribution map.

[0093] The two-dimensional revenue distribution map obtained in step S431 is treated as a continuous image. The circular distribution kernel constructed in step S432 is used as the convolution kernel, and a two-dimensional spatial convolution operation is performed. Specifically, on the global work case grid, for each grid node, the center of the circular distribution kernel is placed at that node. The sum of the products of the original revenue values ​​of all grid points within the kernel function coverage area and the corresponding position weights of the kernel function is calculated as the new revenue value of that node. This process is repeated for all nodes to generate the enhanced revenue distribution map. Because the circular distribution kernel is negative or zero in the center and positive around the edges, the mathematical effect of convolution is equivalent to replacing the original revenue value of each point with the average revenue of its surrounding circular area minus the revenue of its central area. This operation directly results in the enhanced revenue being significantly suppressed near work case points that have been sampled multiple times due to center suppression; while in positions where sampling points are sparse but located on the outer ring of sampling points, the enhanced revenue is significantly amplified due to the high weights from the ring.

[0094] Step S434: Perform local stability screening on the enhanced return distribution map, remove pseudo-extreme patches where the return value changes drastically with small positional shifts, and retain stable extreme connected regions that are robust to positional perturbations.

[0095] After spatial modulation, the enhanced return distribution map may exhibit many sharp, isolated local peaks. These peaks may be spurious extrema caused by numerical noise or kernel function boundary effects, and are not robust. Local stability screening is performed by evaluating the curvature characteristics of each local extremum point. For each local maximum point on the enhanced return distribution map, its Herze matrix is ​​calculated, and the eigenvalues ​​of the Herze matrix are solved. If the minimum eigenvalue is negative and its absolute value is large, it indicates that the extremum point drops sharply in a certain direction, is very sharp, and the return jumps drastically with a slight positional shift, belonging to an unstable structure, and is thus discarded. Only those extrema points where all eigenvalues ​​of the Herze matrix are negative, and the ratio of the maximum to the minimum eigenvalue is less than a stability threshold, are retained. These points and their surrounding contour regions constitute stable extremum connected regions. The stability threshold can be set according to the smoothness requirements of the performance proxy surface, thereby ensuring that the retained regions are all robust to positional perturbations and have broad and gentle return plateau regions.

[0096] Step S435: Overlay and compare the retained stable extreme value connected regions with the feasible region risk penalty field, remove the sub-regions in the stable extreme value connected regions that fall into the strong response region of the feasible region risk penalty field, and obtain the candidate information collection benefit region.

[0097] Stability screening ensures the geometric robustness of the gain region but does not consider its physical feasibility risk. The stable extreme connected regions retained in step S434 are superimposed and compared with the feasible region risk penalty field generated in step S410. For each stable extreme connected region, the feasible region risk penalty field value at the grid points it covers is checked. A risk tolerance threshold is set; if the penalty field value at any grid point within the connected region exceeds the risk tolerance threshold, it indicates that this region is too close to the thermodynamic constraint boundary and poses a testing risk. These sub-regions with excessive risk are directly removed from the stable extreme connected regions. The removal operation divides or shrinks the original connected region. The regions remaining after removal, which simultaneously satisfy the location robustness and physical feasibility risk constraints, are formally identified as candidate information collection gain regions.

[0098] Step S436: Use the enhanced revenue signal value of each point within the candidate information collection revenue area as the corresponding point value of the information collection revenue distribution, and assign zero revenue value to points outside the candidate information collection revenue area to obtain the sparse information collection revenue distribution.

[0099] The process of generating the final information gathering benefit distribution focuses entirely on candidate regions. A benefit distribution map of the same size as the global grid is created and initialized to all zeros. Then, all grid nodes designated as candidate information gathering benefit regions are traversed, and the enhanced benefit signal values ​​calculated in step S433 are copied verbatim to their corresponding positions in the new distribution map. For all points falling outside the candidate information gathering benefit regions, their values ​​remain zero. This results in a sparse information gathering benefit distribution, which has positive values ​​only in a few rigorously selected candidate regions, while being zero in the vast explored region and boundary risk region. This significantly reduces the search space for subsequent extreme value searches and improves decision-making efficiency.

[0100] Step S440: Convert the information collection benefit distribution and the feasible domain risk penalty field into dimensionless information collection utility distribution and risk penalty utility distribution respectively through a preset utility scale mapping function. Then, combine the two point by point, and subtract the corresponding point response of the risk penalty utility distribution from the information collection utility distribution point by point to generate a joint evaluation value distribution.

[0101] The information collection benefit distribution and the feasible region risk penalty field have different physical dimensions and numerical scales, making direct combination and comparison impossible. The pre-defined utility scale mapping function is a set of monotonically increasing nonlinear functions, which uniformly map physically meaningful numerical values ​​to a dimensionless utility scale within the 0-1 interval. For example, the utility mapping function for information collection benefits can be a logarithmic flexible maximum function, allowing linear mapping of the benefit signal in the low-value region and logarithmic compression in the high-value region; the utility mapping function for the feasible region risk penalty field can be an exponential function, disproportionately imposing high-weight penalties on medium- and high-risk areas. By inputting the information collection benefit distribution point-by-point into its corresponding mapping function, the information collection utility distribution is obtained; similarly, by inputting the feasible region risk penalty field point-by-point into its corresponding mapping function, the risk penalty utility distribution is obtained. Finally, point-by-point combination is performed, subtracting the corresponding risk penalty utility value from the information collection utility value at each grid point. If the difference is less than 0, it is set to 0. The resulting value constitutes the joint evaluation value distribution. Each value on this joint evaluation value distribution comprehensively reflects the expected information benefits and the necessary risk costs of testing at that point. The larger the value, the more suitable the point is as the next test condition.

[0102] Step S450: On the distribution of joint evaluation values, multiple local search paths are started with multiple random starting points. Each search path moves step by step along the upward direction of the joint evaluation values. The step size is adaptively adjusted according to the local change trend of the evaluation values ​​at the location of the search point until the path converges to the local extreme point.

[0103] In one implementation, step S450 may include the following steps S451-S456: Step S451: Randomly select a set of starting points from the joint evaluation value distribution. The selection probability of the starting points is biased towards the non-zero region of the information collection benefit distribution, and the working condition space distance between any two starting points is kept at no less than the preset decorrelation interval to avoid premature overlap of search paths.

[0104] The choice of starting point directly affects the probability of discovering the global optimum. To ensure that the starting point is likely to fall within a valuable region, an importance sampling strategy is adopted. Using the information acquisition reward distribution as a non-normalized probability quality function, biased random sampling is performed in regions where the reward value is non-zero. Points with higher information reward values ​​have a greater probability of being selected as starting points. Simultaneously, to ensure path diversity, a hard exclusion condition is introduced, stipulating that the Euclidean working space distance between any two selected starting points must remain within a range not less than a preset decorrelation interval, which is set to two to three times the average local correlation length of the performance proxy surface. Through repeated sampling and exclusion checks, a set of starting points that meets the conditions and reaches the preset total number of paths is selected. This effectively avoids multiple search paths starting from the same attraction basin, leading to redundant searches and wasted computational resources.

[0105] Step S452: For each starting point, determine the upward direction of the joint evaluation value of its location, and move to the candidate point along the direction with the current step size. If the joint evaluation value of the candidate point is higher than that of the current point, accept the move and take the candidate point as the new current point; otherwise, reduce the step size and reverse the sign of the upward direction to retry.

[0106] The direction of ascent of the joint evaluation value can be determined by calculating the numerical gradient of the point on the distribution of joint evaluation values. The direction pointed to by the gradient vector is the direction of the fastest ascent of the function value. Moving along this ascent direction from the current point by the current step size, we reach a candidate point. We then query the value of the joint evaluation value distribution at the candidate point. If the value of the candidate point is indeed higher than that of the current point, the move is considered successful, the candidate point is updated to the new current point, and the step size is appropriately increased to accelerate subsequent movements. If the value of the candidate point decreases instead of increasing, it indicates that the current step size has exceeded the ascent confidence region, so the move is rejected, the current point remains unchanged, and the step size is reduced by a decay factor, for example, to half its original value. If the reduced step size is too small, we attempt to reverse the ascent direction to test the waters, in case the current point happens to be on an extremely narrow local ridge. This testing and step size self-adjustment mechanism achieves a reliable and efficient climbing process.

[0107] Step S453: During the movement, continuously monitor the decay rate of the joint evaluation value of the path points. When the decay rate continuously exceeds the preset critical change rate, it is determined that the search path has entered the area dominated by the feasible domain risk penalty field, and the extension of the search path is terminated in advance.

[0108] Although the starting points are all chosen in areas with informational benefits, moving along the upward direction may inadvertently enter regions near the boundary where the joint evaluation value is dominated by the risk penalty term. In such regions, the joint evaluation value will drop sharply due to the sudden increase in the utility of the risk penalty. In each move, not only is the evaluation value of the current point compared to the candidate points checked, but the decay rate of the joint evaluation value over the most recent moves is also calculated, i.e., the magnitude of the decrease in evaluation value per unit movement distance. If this decay rate exceeds the preset critical rate of change in two consecutive moves, this is a strong signal indicating that the search path is veering out of the feasible information region and into a high-penalty region dominated by the feasible domain's risk penalty field. At this point, without waiting for natural convergence, the extension of this search path should be immediately determined and terminated, discarding it to avoid providing a test suggestion on the verge of danger due to forcibly searching the boundary.

[0109] Step S454: After all search paths have converged to their respective local extreme points or have been terminated prematurely, the final working condition coordinates of each path are aggregated into a set of candidate extreme points.

[0110] The status of all parallel-initiated local search paths is summarized and checked. For paths that are not prematurely terminated and converge through normal step-by-step movement, their final stopping points, i.e., the working coordinates of local extrema, are extracted. For paths that are prematurely terminated in the early stages of movement, their stopping points are discarded. The coordinates of the final stopping points of all normally converged paths are gathered together to form a candidate extrema set. Each point in this set represents a locally optimal trade-off solution of the joint evaluation function within a certain feasible region.

[0111] Step S455: Perform nonmaximum suppression on the candidate extreme point set, merge redundant extreme points whose spatial distance is less than the local correlation length of the performance proxy surface, and retain only the independent candidate extreme points that are spatially isolated from each other.

[0112] The candidate extremum set may contain multiple points that are very close to each other, actually located near the bottom of the same broad local extremum basin, representing the same extremum, and thus redundant. The spatial metric used for nonmaximum suppression is the local correlation length of the performance proxy surface, which defines the correlation radius within a test case space. Within this correlation radius, function values ​​are considered highly correlated, and multiple independent extrema should not exist. The suppression process uses a greedy algorithm: first, candidate extremum points are sorted from highest to lowest joint evaluation value; then, the point with the highest evaluation value is selected as the first independent candidate extremum point and retained; next, the remaining candidate points are traversed, and if the test case space distance between a candidate point and any retained independent candidate extremum point is less than the local correlation length, the candidate point is considered redundant and discarded; this step is repeated until all candidate points have been processed. Finally, the retained points are independent candidate extrema points that are sufficiently isolated in the test case space and each represents a different performance region.

[0113] Step S456: Compare the joint evaluation values ​​of all independent candidate extreme points, and determine the independent candidate extreme point with the highest joint evaluation value as the working condition for the next test to be performed.

[0114] After non-maximum suppression, a small number of independent candidate extreme points are obtained, each representing the optimal choice in different regions under the trade-off between exploration and risk. To make a final decision, the values ​​of these points on the joint evaluation value distribution are directly read, and a simple maximum comparison is performed. The independent candidate extreme point with the highest joint evaluation value is selected, and its corresponding continuous operating condition coordinates are taken. These coordinates are then determined by the system as the operating condition for the next test to be executed. This decision reflects the test action plan that maximizes the expected total utility from a global perspective under the current cognitive state, i.e., information gain minus risk cost.

[0115] Step S460: Collect all local extreme points reached by the search paths, select the local extreme point with the highest joint evaluation value, and determine the continuous working condition coordinates corresponding to the extreme point as the working condition for the next test to be performed.

[0116] Step S500: Drive the chiller unit under test to run according to the operating conditions of the next test to be performed and collect the corresponding performance observation data. Update the posterior distribution of the performance proxy surface through the collected performance observation data, and repeat the determination and execution of the operating conditions of the next test to be performed until the total global prediction uncertainty of the performance proxy surface meets the preset termination judgment condition.

[0117] As one implementation method, step S500 employs an incremental propagation mechanism based on local hidden state tracing. The update process specifically includes the following steps S510-S560: Step S510: After acquiring the newly collected performance observation data, the working condition location identifier in the new observation data is matched with the existing induced location set to determine the target induced location set affected by the new observation in the working condition space. The neighbor matching is based on the association metric rule defined by the initial kernel morphology description set, so that the matching range is adapted to the local association length of the current proxy surface.

[0118] The first step in incremental updates is to determine the scope of influence of the new observation data. The newly acquired performance observation data includes the operating condition location identifier of the test point, i.e., its operating condition embedding representation. The existing set of induced locations consists of the induced points retained after pruning the current supporting surrogate surface. During neighbor matching, the correlation metric between the operating condition location identifier of the new observation point and the operating condition coordinates of each existing induced location is calculated. This metric is calculated using the embedding kernel defined by the initial kernel morphological descriptor set. A matching range threshold is set, which is equal to the local correlation length of the current surrogate surface. All existing induced locations whose correlation metric with the new observation point exceeds this threshold are identified as the target induced location set significantly affected by the new observation. The use of the correlation metric rule ensures that the delineation of the influence range is strictly consistent with the kernel function characteristics of the current Gaussian process surrogate model, preventing the influence from being incorrectly propagated to regions that the model considers irrelevant.

[0119] Step S520: Extract the current hidden layer state description of each induced position in the target induced position set, compare the performance observation in the new observation with the local prediction mean corresponding to each current hidden layer state description, and generate a matching residual characterizing the degree to which the new observation is explained by the current hidden layer state. The local prediction mean is provided by the local prediction mean vector in the sufficient statistics of each induced position.

[0120] For each induced location in the target induced location set, a local predicted mean vector is extracted from its stored variational posterior sufficient statistics. Then, using the correlation metric between the induced location and the observation point, the performance prediction for the new observation point is calculated in the hidden state of this induced location through the aforementioned decoding network or linear readout head. This predicted value is compared with the newly acquired true performance observation, and the absolute value of the difference between the two is calculated, i.e., the matching residual. The magnitude of the matching residual directly reflects whether the current hidden state description can well explain the new observation data. If the matching residual is small, it indicates that the information carried by the new observation data is already included within the current uncertainty cognition of the induced location; if the matching residual is large, it indicates that the new data contains "surprising" information that challenges the current hidden state description, and the knowledge of this induced location needs to be updated.

[0121] Step S530: For induced positions where the matching residual exceeds the preset tolerance level, perform hidden layer state reinitialization. The initial value of reinitialization is jointly oriented by the performance observation of the new observation and the direction information of the matching residual, so that the reinitialized hidden layer state description is close to the performance trend implied by the new observation. At the same time, for induced positions where the matching residual does not exceed the preset tolerance level, maintain the current hidden layer state description unchanged.

[0122] In step S530, a tolerance threshold for the matching residual is set. This threshold is determined by the standard deviation of the average observation noise of historical models stored in the cross-model knowledge embedding mapping structure. This noise standard deviation is estimated during the meta-task learning phase by statistically analyzing the standard deviation of the repeated measurement error of performance observations under the same operating conditions in all historical model test sequences. It reflects the level of irreducible background noise jointly determined by sensor accuracy and environmental fluctuations. Each induced position in the target induced position set is evaluated one by one. If the matching residual of an induced position is less than the tolerance threshold, the knowledge of that induced position is considered sufficient, and its current hidden state description is frozen. Conversely, if the matching residual exceeds the tolerance threshold, the hidden state reinitialization is triggered. The reinitialization is not a random restart, but a targeted local correction. The newly initialized variational mean vector is set to shift the original mean vector towards the direction of the new observed performance data by a step size controlled by the size of the matching residuals. The step size is the matching residuals multiplied by a learning rate factor. The variational covariance matrix is ​​then expanded back to an initial prior diagonal matrix to reflect the forgetting of knowledge and the openness to new information. In this way, the newly initialized hidden state description is forced to be closer to the state that can better explain the performance trend revealed by the new observation data.

[0123] Step S540: The re-initialized hidden state description of the induced position is propagated along the diffusion path defined by the correlation metric rule of the initial kernel morphology description set to the neighboring induced positions that have not been re-initialized. The conductivity tensor in the correlation metric rule is used to apply a smooth adjustment to the neighboring hidden state description so that the global hidden state manifold of the affected region is maintained.

[0124] Local reinitialization disrupts the smoothness of the global hidden state manifold in the affected region. To maintain global continuity, changes in the local state need to be smoothly propagated along the prior dependency path specified by the model. The propagation process uses the reinitialized induced position as the wave source, determining the diffusion path weights between it and all unreinitialized neighboring induced positions. These weights are given by the conductance tensor in the association metric rule. For each neighboring induced position, based on its conductance weight with the wave source, its current hidden state description is slightly adjusted towards the reinitialized state direction. The adjustment amount is calculated as follows: the new mean of the neighboring position equals the original mean plus the conductance weight multiplied by the change in the wave source mean, multiplied by a decay coefficient. This smooth adjustment propagates hop-by-hop until the effect decays to a negligible level, ultimately ensuring that the hidden state manifold can undergo local deformation while maintaining global consistency after absorbing new observation information.

[0125] Step S550: Refresh the performance proxy surface using the global hidden layer state manifold after state propagation, recalculate the performance prediction mean of each working node in the affected area, and update the prediction uncertainty of each working node in the affected area based on the trace of the local conductivity tensor.

[0126] After state propagation is complete, a globally hidden state manifold that has been locally updated is obtained. The affected region is the subspace of the working conditions covered by the induced positions that participated in state propagation. At each global grid node within this subspace, the decoding and uncertainty calculation process in step S350 is re-executed. The mean performance prediction of that node is refreshed with the new hidden state interpolation results, and its prediction uncertainty is updated based on the trace of the node's local conductivity tensor. For regions unaffected by state propagation, the mean performance prediction and prediction uncertainty of their nodes remain unchanged. This incremental, localized refresh method ensures that the computational cost of updating the surrogate surface is proportional to the influence range of the new observation, rather than proportional to the global grid size, significantly improving the efficiency of iterative iterations.

[0127] Step S560: Record the propagation path of state propagation that occurred during this incremental update and the state change at each induced position along the path to form an incremental propagation trajectory record, which can be used as a reference for differentiating the decay rate of different regions when calculating the total uncertainty of the subsequent global prediction.

[0128] To more precisely characterize the evolution of the surrogate model's knowledge, the propagation dynamics of each incremental update are recorded. Detailed records are kept of which induced locations acted as wave sources triggering propagation, which specific edge paths affected which neighboring induced locations, and the state changes in each dimension of the hidden state description at the endpoint of each path (i.e., the norm of the difference before and after the update). This information is serialized and stored as an incremental propagation trajectory record. These records are used to analyze the decay rate of uncertainty in different regions when calculating the total global prediction uncertainty: regions with frequent large state changes indicate that the model is still undergoing strong learning and adjustment in that region, and its uncertainty decays rapidly; while quiescent regions with no propagation activity for a long time have entered a bottleneck period of uncertainty decay. By referring to these trajectory records, the determination of termination conditions can provide a more detailed understanding of the micro-dynamics of the entire learning process.

[0129] As one implementation method, in step S500, the loop termination judgment logic may specifically include the following steps S570~S5120: Step S570: After updating the performance proxy surface each time based on the newly collected performance observation data, extract the prediction uncertainty of each point on the global operating condition grid, and perform depth correction on the prediction uncertainty based on the depth of each point inside the thermodynamic constraint boundary of the refrigeration cycle, so that the uncertainty contribution of points adjacent to the constraint boundary is attenuated, and the depth-corrected uncertainty distribution is obtained.

[0130] Directly calculating the total global prediction uncertainty may introduce a bias: a large amount of physically unreliable extrapolated uncertainty near the boundary can unreasonably inflate the total value, causing the algorithm to misjudge that the overall model is still unclear, thus excessively delaying termination. To address this, depth correction is introduced. The shortest signed distance from each valid point (i.e., unmarked invalid points) on the global operating condition grid to the thermodynamic constraint boundary of the refrigeration cycle is calculated. A positive distance indicates the point is inside the boundary, and a negative distance indicates the point is outside. Based on this boundary-inside depth, a confidence decay factor between zero and one is calculated; the shallower the depth, the closer to the boundary, the lower the factor value. This can be achieved using a Sigmund function or piecewise linear interpolation. The original prediction uncertainty of each point is multiplied by its corresponding confidence decay factor to obtain the depth-corrected uncertainty distribution. This process artificially attenuates the prediction uncertainty contribution of points near the boundary, making the termination decision more focused on the level of understanding within the core feasible region.

[0131] Step S580: Integrate the depth correction uncertainty distribution over the entire effective operating space, and use the volume of the effective operating space as a normalization reference to obtain the total global prediction uncertainty of the current iteration step.

[0132] The total global prediction uncertainty is a macroscopic scalar indicator that measures the overall cognitive remaining of the model. The deep correction uncertainty distribution obtained in step S570 is discretely integrated over the entire effective operating space grid, i.e., the corrected uncertainty values ​​of all effective nodes are summed. Simultaneously, the total volume of the effective operating space is calculated, which can be approximated by accumulating the number of effective nodes and multiplying by the unit grid volume. The volume-normalized average value is obtained by dividing the sum by the total volume. This volume-averaged deep correction uncertainty is the total global prediction uncertainty of the current iteration step, reflecting the degree of prediction inaccuracy remaining per unit volume within the entire physically feasible core performance region under the current test.

[0133] Step S590: Extract the total amount of global prediction uncertainty recorded in the iterative process to form an uncertainty decay trajectory. Perform trend decomposition on the uncertainty decay trajectory to separate the base decay trend reflecting the overall learning process and the fluctuation adjustment trend reflecting the local exploration effect.

[0134] As testing iterates, a historical record is maintained, storing the total global prediction uncertainty calculated at each step. The curve connecting these points forms the uncertainty decay trajectory. To more clearly determine whether the learning process has reached a bottleneck, trend decomposition is performed on this trajectory curve. A smooth, monotonically decreasing curve can be extracted using local weighted regression scatter smoothing, serving as the base decay trend reflecting the overall learning process. The sequence remaining after subtracting the base trend from the original trajectory serves as the fluctuation adjustment trend reflecting the local exploration effect. The fluctuation adjustment trend typically shows a sharp drop when exploring entirely new high-performance areas, and a slight rebound when exploring areas with minimal gains. Through this separation, changes in the base learning rate can be analyzed independently.

[0135] Step S5100: Monitor the slope change of the base decay trend. When the slope of the base decay trend is continuously lower than the preset rate of change in multiple consecutive iterations, update the termination tolerance limit included in the termination judgment condition to the preset proportion of the current total global prediction uncertainty, thereby achieving adaptive relaxation of the termination condition.

[0136] The slope of the base decay trend directly reflects the rate at which the global uncertainty of the model is reduced by adding new test points. This slope is continuously monitored. Once it is detected that the slope remains above a preset minimum rate of change threshold for M consecutive iterations (e.g., M is set to 5), with its absolute value close to zero, it indicates that the base learning rate has become extremely slow, the uncertainty reduction from further testing is negligible, and the model's capability has reached the upper limit provided by the data or priors. At this point, continuing to use the previously set fixed, potentially overly strict termination tolerance would lead to a meaningless infinite loop. Therefore, an adaptive relaxation mechanism triggers the termination condition, automatically updating the preset termination tolerance to the current total global prediction uncertainty multiplied by a proportional coefficient slightly greater than one, such as 1.05. This adaptive relaxation allows the algorithm to converge promptly when performance improvement stagnates, rather than terminating indefinitely due to the inability to achieve unrealistically high accuracy.

[0137] Step S5110: Compare the total global prediction uncertainty of the current iteration step with the updated termination tolerance limit. When the total global prediction uncertainty is not higher than the termination tolerance limit, send a termination signal to the test process.

[0138] At the end of each iteration, the total global prediction uncertainty of the current step is numerically compared with the current termination tolerance limit in the system. The termination tolerance limit may be an initial preset value or a value updated by adaptive relaxation in step S5100. Once it is found that the total global prediction uncertainty is less than or equal to the termination tolerance limit, it means that the model has reached a preset sufficient level of mastery over the global performance, and a high-priority termination signal is sent to the control loop of the entire testing process.

[0139] Step S5120: After sending the termination signal, freeze the state of the performance proxy surface and output the mean of the global performance prediction at the freezing moment as the final performance description of the chiller unit under test.

[0140] Upon receiving the termination signal, the test loop is immediately interrupted, and no new test conditions are generated or issued. At this moment, all data representing the latest known performance proxy surface in memory, especially the global hidden layer state manifold and parameters, are set to read-only frozen state, prohibiting any subsequent update operations. The performance prediction mean of each node on the global condition grid at the time of freezing is used as the complete performance spectrum inferred from the intelligent optimization test of the chiller unit under test. This spectrum is then organized into a final standardized performance description data file for storage and display, for example, plotted as a full-condition performance coefficient characteristic curve. This final performance description is the core product output by this method.

[0141] As one implementation method, after freezing the performance proxy surface, the method further includes selectively solidifying the path states related to the current test task in the cross-machine knowledge embedding mapping structure, which may specifically include the following steps S5130~S5180: Step S5130: Extract the global hidden layer state manifold evolution snapshot sequence generated during all iterations of this test task. The evolution snapshot sequence records the topology and hidden layer state description of the global hidden layer state manifold after each update in the iteration order.

[0142] Throughout the entire testing task, the global hidden state manifold generated in each step S340 is sequentially recorded as a snapshot. After the task ends, these snapshots form an evolutionary snapshot sequence arranged in iteration time order. Each snapshot not only contains the hidden state description values ​​for each node in the global domain, but also implicitly records the topological structure of the hidden state manifold at that moment, i.e., the intrinsic network of inter-node connections defined by the diffusion conductivity tensor. Extracting this evolutionary snapshot sequence yields a complete time-slice archive of the learning dynamics of this testing task, providing rich material for extracting stable knowledge subsequently.

[0143] Step S5140: Perform phase space reconstruction processing on the evolution snapshot sequence, calculate the state transition entropy on the reconstructed set of phase points, identify the stable phases whose state transition entropy is lower than the global average level, and extract a representative hidden layer state manifold snapshot from the stable phase.

[0144] To pinpoint the mature and stable phase of the model throughout the dynamic process, the phase space of the evolution snapshot sequence is first reconstructed. Each frame's global manifold is flattened into a high-dimensional vector, and these vectors are mapped to points in the phase space using time-delay embedding. The state transition entropy is calculated for this phase point sequence; this entropy quantifies the amount of information required to transition from one frame's state to the next, i.e., the complexity of the model's evolutionary dynamics. The average state transition entropy of the entire sequence is calculated. Then, searching backward from the end of the sequence, a continuous, longest interval is located where the local average state transition entropy is consistently lower than the global average. This indicates that the model's structural changes are minimal at this stage, and it has entered a stable convergence phase. The most central hidden layer state manifold snapshot from this stable phase is selected as a representative example of the mature experience gained in this task.

[0145] Step S5150: Reverse map the captured hidden layer state manifold snapshot back to the output space of the adaptive kernel morphological descriptor fusion path, restore it to the set of stable kernel morphological descriptors that dominate the performance proxy surface morphology in this stable phase, and record its corresponding test task context.

[0146] The snapshot of the global hidden state manifold with the highest stability, extracted in step S5140, is reconstructed into a corresponding set of kernel morphological descriptors using an inverse mapping algorithm. This inverse process is essentially an approximate inverse of morphological compression mapping. Using this manifold snapshot as the target, and fixing the other parts of the Gaussian process surrogate model, a set of kernel morphological descriptors is solved to minimize the Frobenius norm difference between the global hidden state manifold generated by the construction process in step S300 and the target manifold snapshot. This inverse problem can be solved using gradient-based optimization methods. The converged set of kernel morphological descriptors is the stable set of kernel morphological descriptors, representing the kernel morphology that, after sufficient learning, is most effective in characterizing the global performance of this model in this test task. Simultaneously, the context information of this test task is recorded, including a summary of the inherent attribute information of the tested model and the core test condition range.

[0147] Step S5160: Compare the stable kernel morphology description set with the historical kernel morphology description prototypes currently stored in the cross-model knowledge embedding mapping structure to separate the overlapping feature channels that are compatible with the historical prototypes and the difference feature channels that conflict with the historical prototypes.

[0148] The cross-model knowledge embedding mapping structure maintains a continuously growing library containing historical model kernel morphology description prototypes extracted and retained from previous complete test tasks. The stable kernel morphology description set obtained in the current task is compared channel-by-channel with each historical prototype in the library. Similarity calculation is based on the Pearson correlation coefficient between activation modes of each channel. By aggregating the comparison results with all historical prototypes, for each feature channel of the stable kernel morphology description set, if it consistently maintains a high correlation coefficient with at least one corresponding channel of a historical prototype, exceeding an overlap threshold, then the channel is classified as an overlapping feature channel, meaning the knowledge learned in the current task is compatible with historical priors. Conversely, if the correlation coefficient of a channel with all corresponding channels of historical prototypes consistently falls below a difference threshold, it is classified as a differing feature channel, indicating that the channel captures new knowledge specific to the current task that conflicts with historical experience. Channels falling between these two thresholds are labeled as neutral feature channels.

[0149] Step S5170: Apply persistence enhancement processing to the modulation descriptors associated with the overlapping feature channels in the unit structure coding path and the operating condition coding path to reduce their natural decay rate in subsequent tasks.

[0150] Overlapping feature channels represent general knowledge that has been repeatedly validated across tasks and should be more firmly remembered by the model. Through channel function tracing, we trace back to the unit structure coding path and operating condition coding path to identify modulation descriptors with strong mapping relationships to these overlapping feature channels. We then apply persistent reinforcement to these located modulation descriptors. Specifically, in the optimizer internal states associated with these modulation descriptors, we increase the cumulative exponential moving average of their historical gradient squares. This is equivalent to enhancing the consolidation coefficients of these parameters in subsequent meta-task learning or resilient weight consolidation algorithms for continuous learning. This significantly limits the allowed update step size and deviation magnitude of these parameters when the model faces new test tasks in the future, effectively reducing the rate of natural forgetting.

[0151] Step S5180: Apply initial response soft masking to the modulation descriptors associated with the difference feature channels in the unit structure coding path and the operating condition coding path to reduce their influence in the cold start phase of new tasks, and persistently store the processed path state.

[0152] In one implementation, step S5180 may specifically include the following steps S5181 to S5185: Step S5181: Calculate the activity level indicator of each feature channel in the stable kernel morphology descriptor set, and extract the activity level indicator of the corresponding channel in the historical model kernel morphology descriptor prototype. Compare the two channels one by one to generate a channel overlap index sequence.

[0153] The activity level indicator for each feature channel in the stable kernel morphology descriptor set can be defined by calculating the variance of the activation values ​​of that channel across all time steps or pseudo-time steps. Similarly, for each historical model kernel morphology descriptor prototype, the activity level indicator for each corresponding channel is also calculated. For each feature channel i, the current activity level indicator is compared one by one with the corresponding channel activity level indicators of all historical model prototypes in the library, by calculating the cosine similarity between each pair of channels. The maximum similarity between the channel and all historical prototypes is taken as the overlap index of that channel. After traversing all channels, a channel overlap index sequence is generated, where each element quantifies the prevalence of the corresponding channel's features being reproduced in historical experience.

[0154] Step S5182: Mark channels in the channel overlap index sequence that exceed the first judgment limit as overlapping feature channels, mark channels that are below the second judgment limit as differential feature channels, and mark channels that are between the first and second judgment limits as neutral feature channels.

[0155] A first threshold, such as 0.8, is set as the first judgment threshold, and a second threshold, such as 0.3, is set as the second. The channel overlap index sequence is then evaluated element-by-element: channels with an index value greater than or equal to 0.8 are marked as overlapping characteristic channels; channels with an index value less than or equal to 0.3 are marked as differential characteristic channels; and channels with an index value between 0.3 and 0.8 are marked as neutral characteristic channels. This simple and clear process categorizes all channels into three classes based on their novelty and generality, providing a precise index for differentiated post-processing.

[0156] Step S5183: For channels marked as overlapping feature channels, trace the modulation descriptors that have a mapping relationship with the channel in the unit structure coding path and the operating condition coding path, and slow down the forgetting rate of these modulation descriptors so that their state can be maintained for a longer time in training across different tasks.

[0157] By using a pre-established index table, it is possible to reverse-engineer which neurons or columns of weight matrices in the graph neural network layers and fully connected layers of the unit structure encoding path, as well as in the dilated convolution and attention layers of the operating condition encoding path, have a strong mapping relationship with the outputs of these labeled overlapping feature channels. These located weights and biases are the associated modulation descriptors. Slowing down their forgetting rate is specifically implemented by modifying the optimizer state of these parameters, increasing the momentum term coefficient decay factor, or adding a secondary penalty term to the loss function for changes in these parameters, with a relatively large penalty coefficient. This measure makes it difficult for the state of these modulation descriptors to shift significantly during training for any subsequent new task, thus enabling this general knowledge to be maintained for a long time across different test tasks.

[0158] Step S5184: For channels marked as differential feature channels, trace the modulation descriptors that have a mapping relationship with the channel in the unit structure coding path and the operating condition coding path, and apply a response attenuation adjustment to these modulation descriptors so that their initial activation level is compressed when facing new tasks.

[0159] Similarly, modulation descriptors strongly correlated with the differential feature channels are traced back using an index table. Response attenuation can be applied by permanently cascading a learnable scaling-gated unit at the output of the modulation descriptor. The bias term of this unit is initialized to a large negative value, such as -2.0, so that even if the modulation descriptor's output is large, its effective output is compressed to a small interval close to zero after passing through the scaling-gated activation function. Simultaneously, the learning rate of the scaling-gated unit's weight parameters is set to a slightly larger value during training for subsequent new tasks, allowing the model to quickly remove this suppression when it discovers the usefulness of the differential knowledge. This soft-masking mechanism effectively reduces the interference that specialized knowledge might cause during the cold start phase of new tasks without discarding any knowledge.

[0160] Step S5185: For the modulation descriptors associated with the channels marked as neutral feature channels, retain their existing update rules and response characteristics without additional intervention. After completing the category processing of all channels, write the adjusted forgetting rhythm and initial activation level information into the state storage area of ​​the corresponding path in the cross-model knowledge embedding mapping structure.

[0161] For channels classified as neutral feature channels whose generality or specificity has not yet been clearly determined, a conservative strategy is adopted. All associated modulation descriptors maintain standard optimizer configurations and activation parameters, without any additional intervention to accelerate forgetting or compress activation, allowing them to naturally update and evolve in new task data. After completing the differential processing of all three types of channel-associated modulation descriptors, all modified forgetting rhythm parameters, such as momentum decay coefficients, loss weight penalty coefficients, and initial activation information, such as scaling gate bias values ​​and learning rate factors, are uniformly serialized and mapped to the state storage area of ​​the corresponding path in the cross-model knowledge embedding mapping structure. This storage area is from the same source as the mapping structure retrieved in step S100, meaning that the next time this method is called to perform performance testing on another new model of chiller unit, the retrieved structure will be the enhanced cross-model knowledge embedding mapping structure, which has been selectively solidified and become smarter, thus achieving continuous accumulation and evolution of knowledge throughout its entire lifecycle.

[0162] Figure 3 This is a hardware entity diagram of a working condition self-optimization decision-making system provided in an embodiment of the present invention, as shown below. Figure 3 As shown, the hardware entity of the self-optimization decision system 1000 includes a processor 1001 and a memory 1002, wherein the memory 1002 stores a computer program that can run on the processor 1001, and the processor 1001 executes the program to implement the steps in the method of any of the above embodiments.

[0163] The memory 1002 stores computer programs that can run on the processor. The memory 1002 is configured to store instructions and applications that can be executed by the processor 1001. It can also cache data to be processed or already processed (e.g., image data, audio data, voice communication data, and video communication data) of the processor 1001 and the various modules in the self-optimization decision system 1000. It can be implemented by flash memory or random access memory (RAM).

[0164] When the processor 1001 executes the program, it implements any of the steps of the chiller performance testing condition self-optimization method based on cross-model knowledge embedding mentioned above. The processor 1001 typically controls the overall operation of the condition self-optimization decision system 1000.

[0165] The above description is merely an embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A self-optimization method for chiller performance testing conditions based on cross-model knowledge embedding, characterized in that, include: The pre-constructed cross-aircraft model knowledge embedding mapping structure is retrieved. The cross-aircraft model knowledge embedding mapping structure is obtained by training through meta-task learning of historical multi-aircraft model test sequences. It includes an operating condition encoding path for mapping the operating condition sequence to the operating condition embedding representation, a unit structure encoding path for mapping the unit's inherent attribute information to the structure embedding representation, and an adaptive kernel morphology descriptor fusion path for fusing the operating condition embedding representation and the structure embedding representation and outputting a set of kernel morphology descriptors. The initial performance observations of the chiller unit under test are obtained under the initial test conditions. The inherent attribute information of the chiller unit under test is retrieved. The inherent attribute information is used as input to generate a structural embedding representation through the unit structure coding path. The initial test conditions are used as input to generate a condition embedding representation through the condition coding path. The structural embedding representation and the condition embedding representation are fused across domains through the adaptive kernel morphology descriptor fusion path. The initial kernel morphology descriptor set is output, and the initial performance observations are used as the initial observation records. A Gaussian process surrogate model is constructed based on the initial set of kernel morphology descriptors. The posterior distribution of the Gaussian process surrogate model is updated by the initial observation records to generate a performance surrogate surface that characterizes the mean and uncertainty of the global performance prediction of the chiller unit under test. Using the performance proxy surface and the preset thermodynamic constraint boundary of the refrigeration cycle, the joint evaluation function consisting of the information acquisition benefit term and the feasible domain risk penalty term is subjected to extreme value search processing to determine the operating conditions of the next test to be performed. The chiller unit under test is driven to run according to the operating conditions of the next test to be performed and the corresponding performance observation data is collected. The posterior distribution of the performance proxy surface is updated by the collected performance observation data, and the determination and execution of the operating conditions of the next test to be performed are repeated until the total global prediction uncertainty of the performance proxy surface meets the preset termination judgment condition.

2. The method according to claim 1, characterized in that, The method involves cross-domain fusion of the structural embedding representation and the operating condition embedding representation through the adaptive kernel morphology descriptor fusion path, outputting an initial kernel morphology descriptor set, and using the initial performance observations as initial observation records, including: The channel attribution and allocation process is performed on the structure embedding representation. The preset refrigeration cycle configuration template set in the cross-model knowledge embedding mapping structure is called. According to the response tendency of each channel of the structure embedding representation to the refrigeration cycle configuration template, the channels are respectively assigned to the structure template channel set and the individual specific channel set. The structure template channel set is combined into a structure template representation, and the individual specific channel set is combined into an individual specific representation. The structure template representation carries the invariant components bound to the refrigeration cycle topology, and the individual specific representation carries the model difference components. The working condition embedded representation is subjected to time-series component separation processing. The main trend representation of the working condition with smooth change is extracted along the working condition sorting direction implied in the working condition embedded representation. The remaining residual after extraction is used as the working condition perturbation representation. The main trend representation of the working condition retains the trend of performance with smooth evolution of working condition. The working condition perturbation representation retains local non-stationary fluctuations. Using the structural template representation as a conditional modulation signal, the modulation descriptor set associated with the structural template representation in the cross-model knowledge embedding mapping structure is retrieved, and the main trend representation of the operating condition is re-modulated channel by channel through the modulation descriptor set, so that the evolution form of the main trend representation of the operating condition is compatible with the cooling cycle configuration, and a steady-state core basis is obtained. Cross-interaction processing is performed on the individual-specific representation and the operating condition disturbance representation. Through the pre-set interaction relationship template in the cross-model knowledge embedding mapping structure, signal exchange and component fusion are performed on the individual-specific representation and the operating condition disturbance representation at corresponding positions to generate disturbance modulation increment that captures the response mode of the operating condition disturbance to the model-specific behavior. The steady-state kernel substrate and the perturbation modulation increment are concatenated along the channel dimension, and the concatenation result is input into the morphological compression mapping. After channel filtering and merging operations, an initial kernel morphological description set with the same dimension as the preset kernel morphological description set is obtained. The transformation rule of the morphological compression mapping is determined by the kernel morphological constraints stored in the cross-model knowledge embedding mapping structure. The initial performance observations, the condition embedding representation corresponding to the initial test condition, and the initial kernel morphology description set are associated and written into the observation cache to obtain the initial observation record containing the mapping relationship between the condition location identifier and the performance observations.

3. The method according to claim 2, characterized in that, The process of performing cross-interaction processing on the individual-specific representation and the operating condition disturbance representation to generate a disturbance modulation increment includes: The individual-specific representation is divided into several disturbance response heads along the channel direction. Before the division, the number of disturbance response heads is determined based on the historical disturbance response diversity indication recorded in the cross-aircraft knowledge embedding mapping structure. Each disturbance response head independently carries a potential disturbance response mode in a type-specific behavior, resulting in a set of disturbance response heads. The operating condition disturbance representation is mapped to generate a disturbance query vector with the same number of disturbance response headers. Each disturbance query vector carries a query description for pattern matching with the corresponding disturbance response header. For each disturbance response head and its corresponding disturbance query vector, the response modes of each channel in the disturbance response head are compared with the disturbance query vector position by position for mode compatibility. Based on the comparison results, the response sensitivity distribution of the disturbance response head to each part of the operating condition disturbance representation is generated. The operating condition disturbance representation is selectively aggregated using each group of response-sensitive distributions. For each disturbance response head, disturbance fragments compatible with the response mode of that head are extracted, and all disturbance fragments are arranged into a conditional disturbance array according to the disturbance response head. The conditional perturbation array and the individual-specific representation are fused position by position. During the fusion process, the positions in the individual-specific representation that are not activated by the conditional perturbation array are subject to biased suppression, while the positions in the individual-specific representation that are activated by the conditional perturbation array are enhanced and preserved, thus obtaining the fused perturbation representation. A time-consistent projection transformation is applied to the fused perturbation representation to map the fused perturbation representation to a time-series arrangement space shared with the steady-state kernel basis, and the perturbation modulation increment is output.

4. The method according to claim 2, characterized in that, The step of concatenating the steady-state kernel basis and the perturbation modulation increment along the channel dimension, and inputting the concatenation result into a morphological compression mapping to generate an initial set of kernel morphological descriptors includes: A numerical range equalization operation is applied to the steady-state kernel basis to adjust the numerical distribution of each channel so that it falls into a consistent distribution range, thereby obtaining an equalized steady-state basis. The distribution range is determined by the steady-state basis distribution specification preset in the cross-model knowledge embedding mapping structure. The perturbation modulation increment is subjected to boundary constraint processing. Based on the perturbation tolerance boundary stored in the embedded mapping structure according to the cross-model knowledge, the response amplitude exceeding the boundary in the perturbation modulation increment is limited to the boundary edge to obtain the corrected perturbation increment. The equilibrium steady-state basis and the correction perturbation increment are connected along the channel direction, and a transition bridging channel is inserted at the connection. The value of the transition bridging channel is jointly determined by the equilibrium steady-state basis and the adjacent channel of the correction perturbation increment, which is used to mitigate the information gap between the two parts, and a joint representation tensor is obtained. The input form of the joint representation tensor is compressed and mapped, and the number of channels is reduced step by step through multiple compression units. Each compression unit contains staggered one-dimensional channel convolution and pointwise nonlinear mapping, and the distribution pattern of active channels is recorded after each compression. At the output of each compression unit, a symmetry adjustment is applied to the current compression representation, and the interrelation structure between each channel is adjusted so that the correlation morphology matrix that can be constructed based on the compression representation falls into the symmetric positive definite domain, thus obtaining a symmetric compression representation. The final first-level symmetric compression representation is output as the initial kernel morphological descriptor set, while the distribution pattern index of the active channels in the last layer of the morphological compression mapping is saved.

5. The method according to claim 1, characterized in that, The process involves constructing a Gaussian process surrogate model based on the initial set of kernel morphology descriptors, updating the posterior distribution of the Gaussian process surrogate model using the initial observation records, and generating a performance surrogate surface characterizing the global performance prediction mean and prediction uncertainty of the chiller unit under test, including: A set of induced positions is configured in the operating space of the chiller unit to be tested, and a hidden state description is assigned to each induced position. The initial value of the hidden state description is determined by the prior form defined by the initial kernel morphology description set, and the initial association between each induced position is determined by the association measurement rule given by the initial kernel morphology description set. Using the working condition location identifiers and performance observations in the initial observation records as conditions, a response likelihood relationship between the hidden layer state descriptions of each induced location and the observation records is generated, and this response likelihood relationship is combined with the prior association morphology derived from the initial set of kernel morphology descriptions to obtain a structured joint distribution. Under the premise of maintaining the conditional independence structure between the hidden state descriptions of each induced position as defined by the prior association morphology, variational inference is performed on the structured joint distribution, and sufficient statistics of each hidden state description are iteratively adjusted so that the lower bound of variational evidence increases in round by round, resulting in a variational posterior hidden state set. Using the variational posterior hidden state set as the source, and the correlation metric rule implied in the prior correlation form as the diffusion conductance tensor, the manifold diffusion is carried out to the continuous working space, and the hidden state at the discrete induction position is propagated to every node of the global working grid, generating a global hidden state manifold. On the global hidden layer state manifold, the mean performance prediction is calculated node by node, and the prediction uncertainty is calculated based on the local diffusion conductivity tensor trace of the node, resulting in a performance proxy surface indexed by the operating condition coordinates. Traverse the performance proxy surface and mark the node prediction uncertainty corresponding to the outer region of the boundary that intersects with the thermodynamic constraint boundary of the refrigeration cycle as invalid.

6. The method according to claim 5, characterized in that, The variational inference of the structured joint distribution, iteratively adjusting the sufficient statistics describing each hidden state, yields a set of variational posterior hidden states, including: Initialize the variational distribution configuration of the hidden layer state description of each induced position, and take the observation samples that are closest to each induced position in terms of working space distance in the initial observation record as local evidence. Use the local evidence to locally warm up the variational distribution configuration so that the variational distribution is initially biased towards local observation. Calculate the mutual information between the hidden state description of each induced position and the global observation sample set under the current variational distribution. Sort the influence of each induced position according to the mutual information and select the induced positions with the highest mutual information to form the active induced set of the current iteration round. For the induction positions in the active induction set, the gradient signal of the sufficient statistic describing the hidden state with respect to the lower bound of the variational evidence is calculated. For the induction positions in the inactive induction set, the sufficient statistic is kept unchanged to maintain the stability of the global structure. The gradient signal is used to perturb and update the sufficient statistics of the active induced position, and a structure preservation constraint is applied after the update to ensure that the relative distance between the hidden state descriptions of the updated active induced position is consistent with the similarity specified by the prior association morphology. Examine the improvement of the variational evidence lower bound after the update. When the improvement is no longer significant, prune inactive induced positions from subsequent calculations and fix the set of induced positions retained after pruning as the new induced layout. Output the final variational posterior hidden state set, where the sufficient statistics describing each hidden state include the inverse representation of the local predicted mean vector and the local correlation matrix at that induced position.

7. The method according to claim 5, characterized in that, The process of updating the posterior distribution of the performance proxy surface using the collected performance observation data, and repeating the determination and execution of the next test conditions until the total global prediction uncertainty of the performance proxy surface meets the preset termination condition, employs an incremental propagation mechanism based on local hidden layer state tracing, including: After acquiring newly collected performance observation data, the working condition location identifier in the new observation data is matched with the existing set of induced locations to determine the set of target induced locations in the working condition space affected by the new observation. The proximity matching is based on the association metric rule defined by the initial kernel morphology descriptor set, so that the matching range is adapted to the local association length of the current proxy surface. Extract the current hidden layer state description of each induced position in the target induced position set, compare the performance observation in the new observation with the local prediction mean corresponding to each current hidden layer state description, and generate a matching residual characterizing the degree to which the new observation is explained by the current hidden layer state. The local prediction mean is provided by the local prediction mean vector in the sufficient statistics of each induced position. For induced positions where the matching residual exceeds the preset tolerance level, the hidden layer state is reinitialized. The initial value of the reinitialization is jointly oriented by the performance observation of the new observation and the direction information of the matching residual, so that the reinitialized hidden layer state description is close to the performance trend implied by the new observation. At the same time, the current hidden layer state description is maintained unchanged for induced positions where the matching residual does not exceed the preset tolerance level. The re-initialized hidden state description of the induced position is propagated to the neighboring induced position that has not been re-initialized along the diffusion path defined by the correlation metric rule of the initial kernel morphology description set. The conductivity tensor in the correlation metric rule is used to apply a smooth adjustment to the neighboring hidden state description so that the global hidden state manifold of the affected region is maintained. The performance proxy surface is refreshed using the global hidden layer state manifold after state propagation, the average performance prediction of each working node in the affected area is recalculated, and the prediction uncertainty of each working node in the affected area is updated based on the trace of the local conductivity tensor. Record the propagation path of state propagation that occurred during this incremental update and the state change at each induced location along the path to form an incremental propagation trajectory record, which can be used as a reference for different regions to differentiate the decay rate when calculating the total uncertainty of the subsequent global prediction.

8. The method according to claim 1, characterized in that, The process involves using the performance proxy surface, combined with a preset thermodynamic constraint boundary of the refrigeration cycle, to perform extreme value search processing on the joint evaluation function consisting of the information acquisition benefit term and the feasible region risk penalty term, to determine the operating conditions for the next test to be performed, including: The thermodynamic constraint boundary of the refrigeration cycle is discretized into a set of boundary points. A repulsive response core is configured for each boundary point. The intensity of each repulsive response core decreases with the distance from the corresponding boundary point. The superposition of all repulsive response cores forms a feasible region risk penalty field on the global operating condition grid. The distribution intensity of the feasible region risk penalty field is determined by the proximity to the constraint boundary. On the performance proxy surface, the predicted mean of each point in the global working condition grid is extracted, and it is compared with the currently discovered global optimal predicted mean. The positive difference part is retained, and the positive difference part is combined with the prediction uncertainty of the point to generate the original information gain signal. The original information gain signal is convolved and modulated with a preset gain pattern modulation function to obtain the information acquisition gain distribution. The gain pattern modulation function has the pattern characteristic of gain decay near the current sampled area and gain enhancement in the moderate distance area. The information collection benefit distribution and the feasible domain risk penalty field are respectively converted into a dimensionless information collection utility distribution and a risk penalty utility distribution through a preset utility scale mapping function. Then, the two are combined point by point, and the corresponding point response of the risk penalty utility distribution is subtracted point by point from the information collection utility distribution to generate a joint evaluation value distribution. On the distribution of the joint evaluation value, multiple local search paths are started from multiple random starting points. Each search path moves step by step along the upward direction of the joint evaluation value. The step size is adaptively adjusted according to the local change trend of the evaluation value at the location of the search point until the path converges to the local extreme point. Collect all local extreme points reached by the search paths, select the local extreme point with the highest joint evaluation value, and determine the continuous working condition coordinates corresponding to the extreme point as the working condition for the next test to be performed.

9. The method according to claim 8, characterized in that, The step of convolving and modulating the original information gain signal with a preset gain pattern modulation function to obtain the information acquisition gain distribution includes: The distribution of the original information revenue signal on the global working condition grid is regarded as a two-dimensional revenue distribution map, and its spatial structure is extracted. The revenue pattern modulation function is constructed as a ring-shaped distribution kernel with a central suppression region and an surrounding enhancement band. The radius of the central suppression region of the ring-shaped distribution kernel is determined by the minimum envelope radius of the currently sampled operating point in the operating space, and the radius of the surrounding enhancement band is determined by the average amount of the local correlation length in the performance proxy surface. The two-dimensional revenue distribution map is spatially convolved with the annular distribution kernel to suppress the original information revenue in the densely sampled working point area and amplify the annular area where the sampled working point is sparse and maintains a moderate distance from the sampled area, thus obtaining the enhanced revenue distribution map. The enhanced return distribution map is subjected to local stability screening to remove pseudo-extreme patches whose return values ​​change drastically with small positional shifts, and retain stable extreme connected regions that are robust to positional perturbations. The retained stable extreme value connected regions are superimposed and compared with the feasible region risk penalty field. Sub-regions in the stable extreme value connected regions that fall into the strong response region of the feasible region risk penalty field are removed to obtain the candidate information collection benefit region. The enhanced revenue signal value of each point within the candidate information collection revenue area is used as the corresponding point value of the information collection revenue distribution, and points outside the candidate information collection revenue area are assigned a zero revenue value to obtain the sparse information collection revenue distribution.

10. A condition-based self-optimization decision-making system, comprising a memory and a processor, wherein the memory stores a computer program executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method according to any one of claims 1 to 9.

Images