Cosmetic stability test-oriented accelerated aging data time series trend prediction method

Abstract

This invention relates to a method for predicting the time-series trends of accelerated aging data for cosmetic stability testing, aiming to address the difficulty in quantifying and predicting the degradation mechanisms of complex formulation components under multiple environmental stresses. The core solution includes: constructing an experimental dataset of active ingredient degradation under multiple stress conditions; after standardization and multidimensional tensor encapsulation; applying symbolic regression to automatically generate an interpretable kinetic rule set; defining a component degradation correlation graph; and dynamically incorporating physical regulatory factors such as temperature, pH, and light to achieve a differentiable graph neural network structure with physical priors. The model employs a hierarchical message passing and gating mechanism, integrating physical rule parameters for multi-step prediction of stability indicators, and supports cross-component migration and fine-tuning of new formulations. This solution improves the accuracy and physical interpretability of automatic deduction of the aging process of cosmetic formulations, enables the location and mechanism attribution of key failure components, and significantly improves the efficiency of stability assessment.

Description

Technical Field

[0001] This invention relates to the field of cosmetic stability prediction and physical mechanism modeling technology, and in particular to a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing. Background Technology

[0002] Currently, cosmetic stability prediction is mainly based on accelerated aging experimental data under multi-stress conditions. Mainstream technical solutions fall into two categories: one is purely data-driven time-series prediction, such as deep learning models like LSTM, GRU, and TCN, which directly extrapolate trends from experimental data; the other is mechanism-based kinetic modeling, such as fitting parameters using classical reaction kinetic equations. Some cutting-edge research attempts to combine the two, forming a hybrid model of "physical mechanism + machine learning," aiming to improve the interpretability and generalization of predictions.

[0003] However, existing technologies have significant shortcomings. Regarding hybrid models, physical mechanisms are often embedded through hard constraints or simple superposition, making it difficult to achieve co-evolution with data-driven approaches. This results in complex training, slow convergence, and redundant architecture, with weak adaptability to new formulations and stress scenarios, as well as weak multi-scale attribution capabilities. Pure data-driven methods, on the other hand, neglect physical laws such as reaction order and activation energy, making them prone to overfitting and poor generalization under new conditions. Their "black box" nature also severely restricts their application in practical scenarios such as compliance and formulation optimization. Summary of the Invention

[0004] This application provides a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing, aiming to solve one of the problems or issues of the prior art mentioned in the background section above.

[0005] The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing provided in this application specifically includes: S1: Obtain the time-series concentration data of each active ingredient in the typical cosmetic ingredient library under multiple stress conditions in accelerated aging experiments, and record the corresponding environmental stress parameter labels to form the original multidimensional dataset; S2: Perform symbolic regression on the original multidimensional dataset to generate an interpretable symbolic dynamics rule set; S3: Using the interpretable symbolic dynamics rule set to define the graph topology, generate a component degradation correlation graph; S4: Calculate the dynamic adjustment factor according to the interpretable symbolic dynamic rule set, and assign the dynamic adjustment factor to the directed edge weight of the corresponding type label in the component degradation association graph to generate the physical guidance graph topology. S5: Construct the backbone of the graph neural network based on the physical guidance graph topology, use a hierarchical message passing mechanism to perform sparse multiplication of the adjacency matrix and gated aggregation operations, and convert the symbolic rule parameters into linear combination coefficients to output the component state hidden vector sequence; S6: Input the latent vector sequence of the component states into the bi-objective loss function for backpropagation optimization to generate a graph neural network prediction model; S7: For components not seen in the new formula, reuse the embedding space of existing component nodes, fine-tune the sign rule adaptation coefficients between them and adjacent nodes through a small number of samples, and generate adaptive prediction model instances. S8: Using the aforementioned adaptive prediction model instance, the accelerated aging process of the current cosmetics is simulated, and the predicted values ​​of stability indicators and component-level attribution analysis results for future time steps are output.

[0006] The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing provided in this application has the following beneficial effects: This application achieves high accuracy and interpretability in cosmetic stability prediction through a collaborative modeling mechanism of symbolic regression and graph neural networks. First, symbolic regression is used to automatically mine kinetic rules with clear chemical meaning that conform to Occam's razor from accelerated aging data, and these rules are transformed into topology and weight adjustment factors for the graph network. This allows the model structure itself to carry physical mechanisms, significantly improving the reliability of predictions across different operating conditions. Second, a lightweight hierarchical message-passing architecture guided in real-time by symbolic rule parameters is designed, employing a dual-objective loss function that integrates multi-step prediction errors and rule consistency constraints. This ensures computational efficiency while forcing the internal state of the network to align with external physical laws, achieving accurate and attributable predictions. Finally, for new formulations, model migration can be quickly completed by reusing existing ingredient embeddings and fine-tuning their symbolic rule adaptation coefficients with neighboring nodes. This significantly reduces the need for new data and deployment costs, providing efficient and reliable intelligent support for formulation development and quality control.

[0007] In summary, this approach achieves a technological breakthrough by deeply integrating symbolic regression and graph neural networks to automatically extract interpretable aging rules from data and guide model structure construction. Without increasing model complexity, it simultaneously improves prediction accuracy, attribution ability, and cross-scenario transfer robustness, forming a new paradigm for cosmetic stability assessment that combines scientific rigor with engineering practicality. Attached Figure Description

[0008] Figure 1 This is the main flowchart of a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing; Figure 2This is a sub-flowchart of a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing; Figure 3 This is another sub-flowchart of a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing. Detailed Implementation

[0009] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0010] The following disclosure provides many different embodiments or examples for implementing different structures of the invention. To simplify the disclosure, specific examples of components and arrangements are described below. Of course, these are merely examples and are not intended to limit the invention. Furthermore, reference numerals and / or letters may be repeated in different examples; such repetition is for simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or arrangements discussed.

[0011] like Figure 1 As shown, this application provides a method for predicting the time-series trend of accelerated aging data for cosmetic stability testing, specifically including: S1: Obtain the time-series concentration data of each active ingredient in a typical cosmetic ingredient library under multiple stress conditions during accelerated aging experiments, and record the corresponding environmental stress parameter labels to form a raw multidimensional dataset containing information on ingredient types and degradation pathways. The multiple stress conditions include: heat, light, oxygen, and pH.

[0012] S2: Perform symbolic regression on the original multidimensional dataset to generate an interpretable symbolic kinetic rule set. Specifically, perform symbolic regression on the original multidimensional dataset to automatically search for a family of kinetic expressions that satisfy the least Occam's razor principle, in order to generate an interpretable symbolic kinetic rule set characterizing key reaction orders, activation energy dependencies, and pH-temperature coupling effects.

[0013] S3: Utilize the interpretable symbolic kinetic rule set to define a graph topology and generate a component degradation correlation graph. Specifically, each active ingredient and its key degradation products in the typical cosmetic ingredient library are mapped as graph nodes, and directed edges with reaction type labels are constructed based on the reaction relationships identified in the interpretable symbolic kinetic rule set. The reaction types include acid-catalyzed hydrolysis, photosensitive oxidation, and metal ion coordination acceleration.

[0014] S4: Calculate the dynamic adjustment factor based on the temperature-sensitive term, pH response term, and light intensity coefficient in the interpretable symbolic kinetic rule set, and assign the dynamic adjustment factor to the directed edge weight of the corresponding type label in the component degradation association graph to generate a parameterized differentiable physical guidance graph topology.

[0015] S5: Construct the backbone of a graph neural network based on the parameterized differentiable physical guidance graph topology, and use a hierarchical message passing mechanism to perform sparse multiplication of adjacency matrices and gated aggregation operations. Use the linear combination coefficients generated in real time by the symbolic rule parameters as message function constraints to output the component state hidden vector sequence explicitly guided by the aging mechanism.

[0016] S6: Input the latent vector sequence of the component states into a bi-objective loss function for backpropagation optimization to generate a graph neural network prediction model with physical consistency correction. The bi-objective loss function includes: a main task loss, used to calculate the multi-step prediction error of concentration, color, and viscosity stability indices; and an auxiliary task consistency regularization term, used to calculate the consistency regularization term of the sign rule parameters to force the evolution trajectory of the edge weights to maintain the same direction as the numerical trend of the dynamic expression.

[0017] S7: For components not seen in new formulations, reuse the embedding space of existing component nodes, fine-tune the symbol rule adaptation coefficients between them and adjacent nodes using a small number of samples, and generate adaptive prediction model instances with cross-formula transfer capabilities.

[0018] S8: Using the aforementioned adaptive prediction model instance, the accelerated aging process of the current cosmetic is extrapolated, and the predicted values ​​of stability indicators for future time steps and the results of component-level attribution analysis are output to complete the lightweight graph neural network prediction process guided by symbolic regression based on aging dynamics.

[0019] S1: Obtain the time-series concentration data of each active ingredient in a typical cosmetic ingredient library under multiple stress conditions during accelerated aging experiments, and record the corresponding environmental stress parameter labels to form a raw multidimensional dataset containing information on ingredient types and degradation pathways. The multiple stress conditions include: heat, light, oxygen, and pH. Specifically, they include: S1.1: Standardize the coding of preservatives, antioxidants, emulsifiers and functional ingredients in the typical cosmetic ingredient library to generate an index set of ingredient types with unique identifiers, which will serve as the basic mapping entity for subsequent multi-stress experimental data association.

[0020] The original ingredient list from a typical cosmetic ingredient library, covering four major categories of chemical substances—preservatives, antioxidants, emulsifiers, and active ingredients—is obtained as input for standardized coding processing. Chemical structure analysis is performed on each category of chemical ingredient in the original ingredient list to extract its International Union of Pure and Applied Chemistry (IUPAC) nomenclature, Chemical Abstracts Service (CAS) Registry Number, and Simplified Molecular Linear Input Specification (SMILES) string, constructing an ingredient feature metadata table containing basic chemical attributes. Based on this ingredient feature metadata table, a multi-level classification index system is established, subdividing preservatives into subcategories such as phenoxyethanol and parabens, and antioxidants into subcategories such as vitamin E derivatives and polyphenols. A unique hierarchical classification code is assigned to each subcategory, forming the first-level semantic index.

[0021] For homologous components with the same core framework but different side chains, topological similarity is calculated using molecular structure feature fingerprints. Components with highly similar structures are grouped into the same structural cluster based on a similarity threshold, and a unique structural cluster identifier is generated for each cluster, forming a second-level structural index. Combining the first-level semantic index and the second-level structural index, a globally unique identifier (UID) is generated for each specific component instance in the component database. This identifier is composed of a category code, a structural cluster ID, and a sequence number, ensuring unambiguous mapping of component entities in subsequent multi-stress experimental data association. The generated globally unique identifier (UID) is bound to the corresponding SMILES string, molecular weight, polar surface area, and other physicochemical parameters, encapsulated into a standardized component category index data structure. This data structure supports fast retrieval based on hash tables and association queries based on graph structures.

[0022] S1.2: Based on the component type index set, perform accelerated aging experiments with multiple stress coupling of heat, light, oxygen and pH, and use high performance liquid chromatography to perform high-frequency sampling and detection of the concentration of active ingredients during the reaction process to generate an original concentration time sequence signal that reflects the dynamic changes of chemical degradation.

[0023] The ingredient index set is read to extract the chemical structure identifiers and physicochemical property parameters of the active ingredients to be tested, serving as the sample configuration basis for the multi-stress coupled accelerated aging experiment. Based on the pre-set experimental design matrix, a multi-dimensional stress combination condition is constructed, including temperature gradient, light intensity level, oxygen partial pressure level, and pH range. Each condition parameter is mapped to the control interface of the environmental simulation chamber. Standardized coded samples of typical cosmetic ingredients are placed in the environmental simulation chamber, and the accelerated aging program is initiated according to the set stress combination condition, ensuring that the samples receive uniform heat, light, oxygen, and pH stimulation in a closed or semi-closed system. During the accelerated aging process, a high-frequency sampling time window is set based on the kinetic reaction rate prediction model. A trace amount of aliquots sample is extracted from the reaction system using an autosampler, and the local chemical reaction is immediately terminated using low-temperature quenching technology to lock in the instantaneous concentration state.

[0024] Pretreatment procedures were performed on the collected trace samples, including protein precipitation removal, organic solvent extraction, and microporous membrane filtration, to eliminate matrix interference and prepare a clear solution that met the injection requirements of high-performance liquid chromatography (HPLC). The pretreated sample was injected into the HPLC system, and separation was performed using a reversed-phase C18 column. A UV-Vis detector or mass spectrometer was used to acquire the chromatographic peak area signals of each active ingredient and its key degradation products at specific wavelengths or mass-to-charge ratios. Calibration curves were performed using external or internal standard methods, converting the chromatographic peak area signals into absolute concentration values. A correspondence between timestamps and concentration values ​​was established, forming a discrete sequence of concentration monitoring data points. Real-time temperature, light irradiance, oxygen concentration, and pH values ​​within the environmental simulation chamber were recorded simultaneously at each sampling time, generating an environmental stress parameter log strictly time-aligned with the concentration data.

[0025] By linking data, the discrete concentration monitoring data point sequence is fused with the environmental stress parameter log to construct the original concentration time series signal sequence that reflects the dynamic changes of chemical degradation. This sequence completely preserves the composition evolution trajectory under multi-stress coupling conditions.

[0026] S1.3: Perform outlier removal and missing value imputation preprocessing on the original concentration time series signal sequence, and simultaneously record environmental stress parameters such as temperature setpoint, light intensity, oxygen partial pressure and solution pH during the experiment to generate a cleaned and corrected component concentration time series data stream and environmental stress parameter label group.

[0027] The system receives the raw concentration time-series signal sequence reflecting the dynamic changes in chemical degradation, output from step S1.2, along with the synchronously acquired raw monitoring log data of the experimental environment. The raw concentration time-series signal sequence may contain outliers caused by fluctuations in high-performance liquid chromatography (HPLC) detection or local data gaps due to uneven sampling intervals, requiring cleaning and correction. A sliding window-based statistical outlier detection method is performed on the raw concentration time-series signal sequence, setting the window length to 5 times the current time step. The mean and standard deviation of the data within the window are calculated. Values ​​exceeding the mean plus or minus 3 times the standard deviation are identified as suspected outliers, and their timestamp index positions are marked to form an outlier index set, thus eliminating interference from instrument transient noise on the kinetic trend fitting.

[0028] For the marked outlier index locations, local polynomial regression interpolation is used for repair. Three valid sampling points before and after the outlier are selected as the fitting basis to construct a second-order polynomial model. The theoretical concentration estimate at the outlier time is calculated and replaced with the outlier values ​​in the original sequence, generating a pre-denoised concentration time-series data stream to ensure the smoothness and physical continuity of the data curve. Missing value checks are performed on the pre-denoised concentration time-series data stream to identify data gaps caused by equipment failure or sampling omissions. For missing segments shorter than two sampling intervals, linear interpolation is used to fill them; for missing segments longer than two sampling intervals, cubic spline interpolation is used to smoothly reconstruct the data based on the slope and curvature characteristics of adjacent valid data points, generating a complete and intact component concentration time-series data stream.

[0029] The system synchronously reads raw records of environmental stress parameters during the experiment, including temperature setpoints, light intensity, oxygen partial pressure, and solution pH. It performs timestamp alignment verification on the environmental stress parameters, mapping the discrete environmental parameters to the corresponding time points in the component concentration time-series data stream. For environmental parameters with slight time deviations, nearest neighbor interpolation is used for matching, ensuring that each concentration sampling time has a unique and accurate description of the environmental stress state. The aligned environmental stress parameters undergo dimensional standardization preprocessing to eliminate numerical differences between different physical quantities. Temperature parameters are converted to the Kelvin scale, light intensity is standardized to milliwatts per square centimeter (m² / cm²) of ultraviolet irradiance, oxygen partial pressure is converted to kilopascals (kPa), and pH remains dimensionless. The standardized environmental parameters are encapsulated into multidimensional vectors corresponding to environmental stress parameter tag groups, forming structured environmental stress parameter tag groups. The cleaned and corrected component concentration time-series data stream is then structurally bound to the generated environmental stress parameter tag groups. Using timestamps as unique keys, a one-to-one correspondence between concentration data and environmental stress data is established, generating standardized data units containing complete temporal evolution information and multiple stress boundary conditions, which serve as the input basis for subsequent spatiotemporal alignment operations of multi-source data.

[0030] S1.4: Based on the component type index set and the environmental stress parameter label group, perform multi-source data spatiotemporal alignment operation to reassemble the discrete sampled component concentration time series data stream according to a unified timestamp and environmental conditions to generate a component-stress correlation data matrix with clear stress boundary conditions.

[0031] S1.5: Perform multidimensional tensor encapsulation processing on the component-stress correlation data matrix to integrate component type identifiers, degradation product intermediate state information, concentration evolution trajectory and corresponding multi-stress environment labels into structured data units to generate an original multidimensional dataset containing complete component type and degradation path information.

[0032] S2: Perform symbolic regression on the original multidimensional dataset to generate an interpretable symbolic kinetic rule set. Specifically, perform symbolic regression on the original multidimensional dataset to automatically search for a family of kinetic expressions that satisfy the least Occam's razor principle, in order to generate an interpretable symbolic kinetic rule set characterizing key reaction orders, activation energy dependencies, and pH-temperature coupling effects. This includes: S2.1: Multi-scale normalization and noise suppression processing are performed on the concentration time series data and environmental stress parameter labels in the original multidimensional dataset to eliminate dimensional differences and extract high signal-to-noise ratio component degradation feature sequences as the standardized input benchmark for the symbolic regression engine.

[0033] S2.2: Based on the component degradation feature sequence, initialize a candidate expression library containing basic mathematical operators and chemical kinetic functions, and use a genetic algorithm to perform random population evolution operations to generate an initial kinetic expression candidate set covering different combinations of reaction orders.

[0034] Based on the high signal-to-noise ratio component degradation feature sequences output from step S2.1, a candidate expression library containing basic mathematical operators and chemical kinetic prior functions is constructed. This candidate expression library serves as the search space for symbolic regression, and its element composition directly determines the physical interpretability and completeness of the subsequently discovered kinetic rules. A set of basic mathematical operators is initialized, including addition, subtraction, multiplication, division, natural logarithm, exponential function, and power operation. These operators are used to construct the skeleton structure of the expressions, ensuring that the generated kinetic equations have the ability to describe nonlinear concentration change trends. Simultaneously, a set of chemical kinetic-specific functions is introduced, specifically including Arrhenius temperature dependence terms, Langmuir adsorption isotherm variants, and Michaelis enzyme kinetic approximations, to explicitly embed the thermal activation and surface reaction mechanisms in the cosmetic aging process. A genetic algorithm is used to initialize a random population, with each individual representing a potential kinetic expression tree structure. The population size is set to N, where N ranges from 500 to 1000 to ensure the diversity of the search space. Each leaf node of the expression tree consists of component concentration variables, time variables, and environmental stress parameters (temperature, pH, light intensity), while the internal nodes consist of the aforementioned basic mathematical operators or chemical kinetic functions. By randomly generating the initial tree structure, it is ensured that the population is uniformly distributed in the solution space, avoiding getting trapped in local optima.

[0035] The evolutionary operations of the genetic algorithm include three stages: selection, crossover, and mutation. In the selection stage, a tournament selection strategy is used to randomly select k individuals from the current population, retaining the individual with the smallest fitting error and lowest structural complexity to enter the mating pool. In the crossover stage, two parent expression trees are randomly selected, their subtree segments are swapped, and new offspring individuals are generated, thus combining feature patterns from different response paths.

[0036] During evolution, the fitness function value of each individual is calculated in real time. The fitness function is a weighted sum of the squared predicted residuals and a complexity penalty term. The fitness value Fit_i of the i-th individual is calculated using the following formula: Where M is the total number of time series data points, C pred,i,j C represents the predicted concentration value of the i-th expression at the j-th time step. act,j Let L(E) be the actual measured concentration value at the j-th time step, and λ be the complexity penalty coefficient. i ) represents the i-th expression E i The total number of nodes, Fit i Let be the fitness value of the i-th individual. This formula ensures that, while pursuing high accuracy, a simple dynamic model that conforms to Occam's razor principle is prioritized.

[0037] After multiple generations of evolutionary iterations, the evolutionary process terminates when the preset maximum number of generations or fitness convergence threshold is reached. The top K individuals with the highest fitness in the population are extracted to form a candidate set of initial kinetic expressions. This candidate set covers different order combinations from zero-order reactions to complex multi-level coupled reactions, as well as various coupling forms of temperature, pH, and light intensity.

[0038] S2.3: Perform multidimensional stress factor fit evaluation and complexity penalty calculation on the initial dynamic expression candidate set to screen out the preferred dynamic expression subset with the smallest residual and simplest structure under multiple stress conditions of heat, light, oxygen and pH.

[0039] For each mathematical expression structure in the initial candidate set of kinetic expressions, the sum of squared residuals under multidimensional stress conditions is calculated to quantify the fitting deviation between the model predictions and the measured concentration data. A fitness evaluation function based on the least squares method is constructed, substituting time-series data under multiple stress conditions (heat, light, oxygen, and pH) into the candidate expressions to calculate the sum of squared differences between the predicted and actual concentrations at each time step. The fitting error index of the i-th candidate expression is calculated using the following formula: Among them, E fit Here, N is the fitting error index, and C is the total number of sampling points. pred For candidate expressions under specific stress conditions (temperature T) j Illumination j pH j At the next time t j The predicted concentration, C obs Let time t j The actual observed concentrations were used. The calculation process traversed all stress combination experimental groups to ensure that the generalization ability of the expression under all operating conditions was accurately evaluated.

[0040] A depth-first traversal of the syntax tree structure of candidate expressions is performed to count the number of operator nodes, the frequency of variable occurrences, and the nesting depth, in order to construct a structural penalty term representing model complexity. An Occam's Razor-like complexity penalty coefficient is introduced to prevent symbolic regression from falling into overfitting local optima; that is, expressions with simple structures and the ability to explain data variations are prioritized. The complexity penalty function is defined as follows: Among them, P comp N is the complexity penalty value, λ is the preset penalty weight coefficient, and N is the complexity penalty value. op N represents the total number of arithmetic operators and function nodes in the expression. var D represents the total number of times the independent variable appears. nestThis represents the maximum nesting depth of the syntax tree. This penalty mechanism suppresses the selection of complex expressions containing redundant higher-order terms or oscillatory terms with no physical meaning.

[0041] The fitting error index and complexity penalty value are linearly weighted and fused to generate a comprehensive fitness score, which serves as the sole criterion for selecting the optimal kinetic expression. A lower comprehensive fitness score indicates that the expression maintains low prediction error while possessing a simpler physical structure and better conforming to the essential characteristics of chemical kinetics. The comprehensive fitness score S... total The calculation logic is the sum of the fitting error and the complexity penalty, i.e., S. total = E fit + P comp The initial set of candidate kinetic expressions is globally sorted according to their comprehensive fitness scores from smallest to largest. The bottom 50% of high-scoring, low-quality expressions are removed, retaining only the high-potential subset. Physical consistency checks are performed on the retained subset, examining the monotonicity of stress factors and their compatibility with the chemical degradation mechanism, eliminating mathematical solutions that violate fundamental thermodynamic laws. For example, it verifies whether the temperature term coefficient is positive to conform to the energy barrier characteristics of the Arrhenius equation, and confirms whether the pH term's trend in acidic or alkaline regions is consistent with known hydrolysis reaction types. Expressions with high statistical fit but absurd physical meanings are eliminated, such as non-physical cases where negative activation energy or increased light intensity leads to a decrease in degradation rate.

[0042] S2.4: Based on the preferred subset of kinetic expressions, analyze the coefficients of the temperature-sensitive term, pH response term, and light intensity coupling coefficient, and optimize the values ​​of each coefficient using the nonlinear least squares method to generate a prototype of parameterized kinetic rules with clear physical meaning.

[0043] Structural terms containing temperature, pH, and light variables were extracted from the preferred subset of kinetic expressions. Exponential terms in the form of the Arrhenius equation, polynomial terms in the form of acid-base catalysis, and linear product terms in the form of photon yield were identified as physical carriers for the coefficients to be optimized. For the identified temperature-sensitive terms, a nonlinear mapping model based on the Arrhenius law was constructed to convert the experimentally recorded temperature sequence into an absolute temperature scale and initialize the predicted values ​​of activation energy and pre-exponential factor, establishing a theoretical correlation between temperature and the reaction rate constant. The Levenberg-Marquardt algorithm was used to iteratively optimize the coefficients of the temperature-sensitive terms. The residual function was defined as the sum of squares of the differences between the measured concentration change rate and the theoretically predicted rate. By calculating the Jacobian matrix to approximate the Hessian matrix, the damping factor was dynamically adjusted to balance the convergence characteristics of the gradient descent method and the Gauss-Newton method.

[0044] Where r(T) is the reaction rate at temperature T, A is the pre-exponential factor, and E... a R is the activation energy, R is the ideal gas constant, and T is the absolute temperature.

[0045] For the pH response term, based on the acid catalysis or base catalysis order determined in the preferred expression, a polynomial coupling model of hydrogen ion concentration and reaction rate is constructed. The measured pH value is converted into hydrogen ion activity, and a buffer capacity correction coefficient is introduced to eliminate the nonlinear deviation caused by ion intensity fluctuations.

[0046] Where k(pH) is the rate constant at a specific pH, k H and k OH These are the rate coefficients for acid catalysis and base catalysis, respectively, where n and m are the reaction orders, [H + ] and [OH - [These represent the concentrations of hydrogen ions and hydroxide ions, respectively.]

[0047] The acid / base catalytic coefficients and reaction order parameters are simultaneously optimized using the nonlinear least squares method. Boundary constraints are set to ensure the rationality of the physical meaning, such as limiting the reaction order to non-negative real numbers. Redundant terms with statistical significance below a preset threshold are eliminated using confidence interval tests. For the light intensity coupling coefficient, a linear or saturated kinetic model of light intensity and photodegradation rate is established. Zero-order or first-order photoreaction mechanisms are selected based on the preferred expression type. Ultraviolet irradiance data are normalized, and photon yield parameters are fitted to quantify the proportional relationship between photon energy conversion into chemical bond breaking efficiency.

[0048] Where v is the photodegradation rate, Φ is the quantum yield, I is the light intensity, ε is the molar absorptivity, and C is the component concentration.

[0049] Joint optimization of multi-stress coupling coefficients is performed, constructing a comprehensive loss function that includes cross-terms for temperature, pH, and illumination. A regularization term is introduced to penalize collinearity among coefficients. The contribution weights of each stress factor are decoupled using the alternating direction multiplier method, ensuring the physical consistency of each individual coefficient under independent stress conditions. The optimized coefficients are then backfilled into the preferred kinetic expression skeleton, generating a complete kinetic equation containing specific numerical parameters. Each parameter is assigned a clear physical unit and dimension, forming an interpretable prototype of parameterized kinetic rules.

[0050] S2.5: Perform semantic mapping and logical verification processing on the parameterized dynamic rule prototype, and transform the verified rules into a standardized, interpretable, symbolic dynamic rule set to fully characterize the key reaction order, activation energy dependence and multi-stress coupling effect.

[0051] The system receives a prototype of a parameterized kinetic rule, which includes temperature-sensitive term coefficients, pH-response term coefficients, and light intensity coupling coefficients optimized using nonlinear least squares, along with the corresponding mathematical expression structure. Semantic parsing is performed on the mathematical expressions in the prototype to extract chemical entity identifiers such as reactants, products, catalysts, and environmental stress factors, constructing a component-reaction-stress triplet mapping relationship. Based on the extracted chemical entity identifiers, standard chemical names and CAS numbers are retrieved from a database of typical cosmetic ingredients. Entity alignment is performed to eliminate semantic ambiguity caused by differences in experimental nomenclature, ensuring that the chemical components involved in the rule have unique and standardized identifiers.

[0052] The derived kinetic expressions undergo physical validity checks, including verifying whether the reaction order is a non-negative real number, whether the activation energy value conforms to the physical range of the Arrhenius equation, and whether the sign of the rate constant satisfies the positive constraint of the law of mass action. For expressions involving multi-stress coupling effects, the monotonicity and boundary conditions of temperature, pH, and light intensity variables are verified, eliminating anomalous rule segments exhibiting non-physical oscillations or divergent behavior under extreme stress conditions. The validated expression structure is then transformed into a standardized symbolic tree structure, where the root node represents the reaction rate, and the child nodes correspond to concentration, temperature correction, pH correction, and light correction terms, forming a clearly hierarchical logical expression.

[0053] Each leaf node in the symbolic tree is assigned a clear physical dimension label, such as molar concentration, Kelvin temperature, dimensionless pH value, and watts per square meter of light intensity, ensuring dimensional consistency in subsequent calculations. Based on the principle of least Occam's razor, redundancy is pruned from the symbolic tree structure, removing secondary stress terms with coefficients close to zero and contributions to prediction error below a preset threshold, simplifying rule complexity while preserving core dynamic characteristics. The simplified symbolic tree is serialized into an interpretable symbolic dynamic rule set in JSON format. Each rule contains a unique rule ID, associated component ID, reaction type label, standardized expression string, and parameter dimension definitions.

[0054] like Figure 2 As shown, S3: A graph topology is defined using the interpretable symbolic kinetic rule set to generate a component degradation correlation graph. Specifically, each active ingredient and its key degradation products in the typical cosmetic ingredient library are mapped as graph nodes, and directed edges with reaction type labels are constructed based on the reaction relationships identified in the interpretable symbolic kinetic rule set. The reaction types include acid-catalyzed hydrolysis, photosensitive oxidation, and metal ion coordination acceleration. Specifically, this includes: S3.1: Based on the list of active ingredients and their key degradation products identified in the interpretable symbolic kinetic rules set, a unique identifier mapping process is performed on each chemical component entity to generate a standardized graph node set containing specific substance instances such as nicotinamide and vitamin E acetate.

[0055] S3.2: Utilize the chemical transformation logic between nodes in the standardized graph node set, and perform directed connection relationship determination operations based on the reaction path described by the interpretable symbolic kinetic rule set to generate a potential edge connection matrix representing the flow direction from raw materials to intermediates and then to the final product.

[0056] The system reads the entity identifiers of each chemical component from the standardized graph node set, as well as the list of reactant-product correspondences contained in the interpretable symbolic kinetic rule set, as the basic input data for constructing graph topological connections. It iterates through each kinetic expression in the interpretable symbolic kinetic rule set, parsing the reactant component identifiers on the left and the product component identifiers on the right, establishing a unidirectional mapping logic from the reaction source to the degradation products. For each parsed reactant-product combination, it searches the standardized graph node set for the existence of a corresponding source node and target node. If both exist, the reaction path is deemed physically feasible in the graph structure. Based on the determined physical feasibility, an N×N sparse adjacency matrix is ​​initialized, where N is the total number of nodes in the standardized graph node set, and all matrix elements are initially set to zero to store potential directed connection states.

[0057] For each reaction path deemed feasible, obtain the index position i of the source node in the node set and the index position j of the target node. Set the element value of the i-th row and j-th column of the adjacency matrix to 1, indicating the existence of a potential directed edge from component i to component j. If there are multiple dynamic rules triggered by different stress conditions or different reaction mechanisms between the same pair of nodes, retain a single connectivity marker in the adjacency matrix to avoid redundant storage of repeated edges and ensure the simplicity of the graph structure.

[0058] S3.3: For each valid connection edge in the potential edge connection matrix, extract the reaction mechanism description information corresponding to the interpretable symbolic kinetic rule set, and perform classification and labeling processing for specific reaction types such as acid-catalyzed hydrolysis, photosensitive oxidation, or metal ion coordination acceleration to generate a reaction type attribute set with clear physical semantic labels.

[0059] For each valid directed edge established in the potential edge connection matrix, pointing from the active ingredient node to the degradation product node, the mathematical expression structure and parameter set corresponding to the reaction path are extracted from the interpretable symbolic kinetic rule set generated in step S2.5. The stress-dependent term function form in the mathematical expression is analyzed to identify whether it contains hydrogen ion concentration power terms, light intensity linear coupling terms, or metal ion concentration coordination terms, which serve as the physical basis for reaction mechanism classification. For reaction paths containing hydrogen ion concentration [H+] or negative pH exponent terms, the dominant mechanism is determined to be acid-catalyzed or base-catalyzed hydrolysis. The coefficient sign and order of the pH response term in the expression are extracted. If the coefficient is positive and accompanied by a high-order power of [H+], it is marked as acid-catalyzed hydrolysis; if the coefficient is negative and accompanied by a related term of [OH-], it is marked as base-catalyzed hydrolysis. The classification result is mapped to the standardized physical semantic label "acid-catalyzed hydrolysis," and the corresponding pH sensitivity threshold range is recorded.

[0060] For reaction pathways whose expressions contain a product term of light intensity I or ultraviolet irradiance E_uv, their dominant mechanism is determined to be photosensitive oxidation or photolysis. The expression is examined to check for linear or nonlinear coupling between the quantum yield constant and light intensity, confirming that the reaction rate's dependence on light intensity conforms to photochemical laws. These reaction pathways are labeled with the physical semantic tag "photosensitive oxidation," and the light intensity response coefficient is recorded to distinguish the different kinetic behaviors of thermal oxidation and photoinduced oxidation. For reaction pathways whose expressions contain a coordination term or complexation constant term for a specific metal ion concentration [M^n+], their dominant mechanism is determined to be metal ion coordination-accelerated degradation. The interaction order between the metal ion term and the main reactant in the expression is analyzed to identify whether there is a synergistic catalytic effect. These reaction pathways are labeled with the physical semantic tag "metal ion coordination acceleration," and the coordination equilibrium constant is retained as a baseline parameter for subsequent weight calculations.

[0061] For complex reaction pathways involving multiple stress dependencies, a main effect contribution assessment is performed. The numerical proportion of each stress term under standard accelerated aging conditions is calculated, and the stress term with the largest contribution is selected as the primary reaction type label, while the remaining stress terms are retained as auxiliary adjustment factors in the edge attributes. If the contributions of each stress term are comparable, a composite type label "multi-stress coupled degradation" is assigned, and the weight ratio of each sub-item is stored in the metadata. The above classification and labeling results are integrated into the attribute fields of the potential edge connection matrix to form a reaction type attribute set with clear physical semantic labels. This attribute set not only contains discrete category labels but also embeds the initial values ​​of key kinetic parameters, such as the activation energy estimate, pH sensitivity index, and photonic efficiency coefficient, ensuring that each label has quantifiable physical support.

[0062] S3.4: Combining the standardized graph node set, potential edge connection matrix, and reaction type attribute set, perform graph structure assembly and topology integrity verification operations to eliminate isolated nodes and invalid loops, so as to generate a component degradation correlation graph with strict physical prior constraints and a complete structure.

[0063] The algorithm receives a set of standardized graph nodes, a potential edge connection matrix, and a set of reaction type attributes as input data to construct an initial topological framework for a component degradation correlation graph with physical prior constraints. For each component node in the standardized graph node set, the algorithm retrieves edge records in the potential edge connection matrix that have in-degree or out-degree connections to that node, and counts the number of adjacent edges for each node to identify isolated nodes. An isolated node determination threshold of zero adjacent edges is set. Intermediate nodes with zero adjacent edges and not terminal stable products are marked as invalid isolated nodes and removed from the standardized graph node set to prevent virtual nodes without reaction paths from interfering with subsequent message passing mechanisms. For the remaining valid node subset, the algorithm extracts the start and end node indices of all directed edges based on the potential edge connection matrix, constructing a temporary directed graph adjacency list structure for executing the loop detection algorithm.

[0064] A depth-first search is used to traverse the adjacency list of a temporary directed graph, maintaining a stack of currently visited paths and a global set of visited node markers. During traversal, the system continuously monitors for back edges pointing to nodes already existing in the current path. When a back edge is detected, the closed path formed by it is considered an invalid chemical cycle loop. This is because the degradation process of cosmetic ingredients typically exhibits an irreversible or unidirectional reaction flow thermodynamically, and a reverse cycle violates the principles of mass conservation and entropy increase. The system records the set of edge indices for all identified invalid loops, removes these edges from the potential edge connection matrix, severing connections that lead to logical paradoxes and ensuring the directed acyclic or quasi-unidirectional flow characteristics of the graph structure. For non-source raw material nodes whose in-degree or out-degree becomes zero due to the removal of loop edges, an isolated node removal operation is recursively performed until no isolated nodes or invalid loops remain in the graph, forming a topologically connected skeleton structure. The physical semantic labels in the reaction type attribute set are mapped to the retained valid edge connections. Based on the start and end node indices of the edges, labels such as acid-catalyzed hydrolysis and photosensitive oxidation are bound to specific edge objects, completing attribute injection. The cleaned and standardized graph node set, the denoised valid edge connection matrix, and the edge attributes bound to physical semantic labels are integrated to assemble the final component degradation association graph data structure.

[0065] For example, when processing cosmetic formulation data containing niacinamide, vitamin E acetate, and phenoxyethanol, the normalized graph node set contains 15 nodes, where node N1 represents niacinamide, node N2 represents nicotinic acid, node N3 represents vitamin E acetate, node N4 represents vitamin E, and node N5 represents phenoxyethanol. The potential edge connection matrix shows that N1 points to N2 (hydrolysis), N3 points to N4 (oxidation), and N5 has neither an out-degree nor an in-degree (assuming stability under experimental conditions and no degradation product connections). First, the number of adjacent edges is counted. It is found that node N5 has both an in-degree and an out-degree of 0, and N5 is not the final stable endpoint to be monitored (or is marked as an inert carrier rather than an active ingredient). Therefore, it is determined to be an isolated node and removed from the node set, leaving 14 valid nodes. Next, an adjacency list is constructed, and a depth-first search is performed. It is assumed that there is an incorrect path in the data: N2 points to N6 (hypothetical intermediate), N6 points to N1, forming a closed loop N1->N2->N6->N1. When the DFS algorithm starts from N1 and attempts to access N6 while also trying to access N1, it finds that N1 is already in the current recursion stack, and determines that N6->N1 is an invalid cycle edge. This edge is removed from the connection matrix. N6 is then checked; if N6 exists solely because of this cycle, then N6 also becomes an isolated node and is removed. Finally, valid unidirectional edges such as N1->N2 and N3->N4 are retained. The label "acid-catalyzed hydrolysis" generated in S3.3 is bound to the N1->N2 edge, and "photosensitive oxidation" is bound to the N3->N4 edge. The final generated component degradation correlation graph contains 12 valid nodes and 8 directed edges with explicit physical labels, eliminating 1 isolated node and 1 invalid cycle, ensuring the physical rationality of the graph structure.

[0066] S3.5: Based on the generated component degradation correlation graph, traverse all directed edges with physical semantic labels, solidify the category information in the reaction type attribute set into edge metadata, and output a physical guidance graph topology structure that can be used for subsequent dynamic adjustment factor assignment.

[0067] like Figure 3 As shown, step S4: Calculate the dynamic adjustment factor based on the temperature-sensitive term, pH-responsive term, and light intensity coefficient in the interpretable symbolic kinetic rule set, and assign the dynamic adjustment factor to the directed edge weights of the corresponding type labels in the component degradation correlation graph to generate a parameterized differentiable physical guidance graph topology. Specifically, this includes: S4.1: Based on the temperature-sensitive terms parsed from the interpretable symbolic kinetic rules set, the Arrhenius activation energy parameter and pre-exponential factor are extracted. Combined with the real-time temperature values ​​in the environmental stress parameter labels of the current accelerated aging experiment concentration time series data, an exponential function operation is performed to generate a temperature dynamic adjustment factor sequence characterizing the thermal stress-driven effect.

[0068] The system receives the topological object of the component degradation correlation graph with physical prior constraints generated in step S3, and the interpretable symbolic kinetic rule set output in step S2. A subset of kinetic rules containing temperature-sensitive terms is selected from this subset. This subset explicitly records the Arrhenius activation energy parameter and pre-exponential factor corresponding to each chemical reaction path, serving as the fundamental physical constants for calculating the thermal stress-driven effect. Each directed edge in the component degradation correlation graph is traversed, and the reaction type label and corresponding kinetic rule index bound to the edge metadata are read. The real-time temperature value of the current time step is extracted from the environmental stress parameter label group. The extracted temperature value is converted to absolute temperature scale units to ensure consistency with the gas constant units in the Arrhenius equation, eliminating calculation biases caused by dimensional differences.

[0069] For each directed edge corresponding to a specific chemical reaction, a pre-defined Arrhenius exponent calculation module is invoked, and the activation energy and pre-exponential factor parameters specific to that reaction are substituted. The temperature dynamic adjustment factor is calculated using the following formula: Where, k T Here, A is the temperature dynamic adjustment factor, and E is the pre-exponential factor. a Let R be the Arrhenius activation energy, R be the ideal gas constant, and T be the real-time absolute temperature.

[0070] Numerical stability checks are performed on the calculated temperature dynamic adjustment factor. Overflow is prevented by limiting the maximum value of the exponential term, and minimum values ​​are truncated to maintain the effectiveness of gradient propagation. The validated temperature dynamic adjustment factor is serialized into a one-dimensional tensor aligned with the time step of the graph nodes, and labeled as the thermal stress response weight component of the directed edge at the current moment.

[0071] S4.2: Using the pH response term function structure identified in the interpretable symbolic kinetic rule set, read the real-time acid and alkalinity values ​​in the environmental stress parameter labels, substitute them into the preset acid catalysis or base catalysis rate equation model for nonlinear mapping calculation, so as to generate a pH dynamic regulation factor sequence characterizing the chemical environment driving effect.

[0072] The real-time pH values ​​recorded in the environmental stress parameter label are read as the raw input variable for calculating the pH dynamic adjustment factor. This value is directly derived from the data stream collected by the online pH sensor during the accelerated aging experiment and cleaned and corrected in step S1.3, ensuring consistency between the physical state and the chemical environment. pH response terms for specific reaction types are extracted from the interpretable symbolic kinetic rule set. For directed edges labeled "acid-catalyzed hydrolysis," an exponential dependence model for hydrogen ion concentration is identified; for directed edges labeled "base-catalyzed hydrolysis," a linear or nonlinear dependence model for hydroxide ion concentration is identified; for components involving zwitterionic equilibrium, a prototype of the generalized acid-base catalytic rate equation, including isoelectric point shift correction, is extracted. The real-time pH values ​​are converted into hydrogen ion activity or hydroxide ion activity. Based on the principle of solution chemical equilibrium, the hydrogen ion concentration is calculated using the following formula: in, This is a real-time pH value. This corresponds to the molar concentration of hydrogen ions. If the reaction mechanism involves base catalysis, then further calculations using the ion product constant of water are needed. Calculate the hydroxide ion concentration: Substituting the pre-defined rate equations for acid or base catalysis into the model, nonlinear mapping calculations are performed. For acid-catalyzed reaction pathways, the mapping function is constructed using the generalized Brønsted catalysis law: in, This is the coefficient of the acid catalytic rate constant. The reaction orders for acid-catalyzed reactions are determined by the symbolic regression engine in step S2. A similar form is used for base-catalyzed reaction pathways: in, This is the coefficient of the base catalytic rate constant. denoted as the order of the base-catalyzed reaction.

[0073] The calculated raw catalytic rate values ​​were normalized to generate a dimensionless sequence of pH dynamic adjustment factors. The catalytic rate under standard reference conditions (e.g., pH=7.0 or the initial pH of the formulation) was selected as the baseline value. The relative adjustment factor is calculated using the following formula: ,in, This is the calculated value of the original catalytic rate under the current pH conditions. The specific value is obtained by converting the real-time monitored pH value into hydrogen ion or hydroxide ion concentration and substituting it into the equation. or: The calculations show that this treatment eliminates dimensional differences between different reaction types, ensures that the adjustment factor varies within the range of 0 to positive infinity, and has a value of 1 near the neutral reference point, facilitating subsequent multiplicative fusion with temperature and light factors.

[0074] The generated pH dynamic regulation factor sequence is spatiotemporally aligned with the directed edges in the component degradation correlation graph. Based on the timestamp index, the relative regulation factor calculated at each time step is assigned to the edge weight attribute field of the corresponding reaction type, forming a pH-driven weight vector that evolves over time.

[0075] Through the aforementioned nonlinear mapping and normalization process, the environmental stress parameters from the previous step are transformed into a sequence of pH dynamic regulation factors characterizing the chemical environment-driven effects. This achieves the expected technical effect of microscopically transforming macroscopic pH measurements into chemical reaction rate modulation coefficients, providing precise chemical kinetic constraints for the subsequent construction of parameterized differentiable physical guidance graph topologies.

[0076] S4.3: Based on the light intensity coefficient defined in the interpretable symbolic dynamics rule set, obtain the real-time ultraviolet irradiance data in the environmental stress parameter label, perform the photon yield and light intensity product operation to generate a sequence of light dynamic adjustment factors characterizing the photo-oxidation driving effect.

[0077] Irradiance coefficients associated with photosensitive oxidation reaction types are extracted from the set of interpretable symbolic kinetic rules. This relationship defines the nonlinear mapping function structure between ultraviolet irradiance and photochemical reaction rate. Real-time ultraviolet irradiance values ​​at the corresponding time step are read from the environmental stress parameter tag group. These values ​​characterize the energy flux density of the light source within a specific wavelength range in the accelerated aging experimental chamber. The read real-time ultraviolet irradiance data undergoes dimensional consistency verification to ensure that its units are consistent with the light intensity reference units defined in the kinetic expression generated by symbolic regression, eliminating numerical deviations caused by sensor calibration differences. Based on the photon yield parameter in the irradiance coefficient, the electronic transition efficiency constant of the current active ingredient under specific wavelength excitation is obtained. This constant reflects the probability of chemical bond breaking or free radical generation after molecule absorption of photons. The product of photon yield and corrected real-time ultraviolet irradiance is performed to construct a photon flux driving factor characterizing the effective photo-oxidation reaction initiation per unit time. The irradiance dynamic adjustment factor is calculated using the following formula: in, Let be the sequence element of the illumination dynamic adjustment factor at time t. The photon yield coefficients are obtained analytically from the set of symbolic dynamical rules. The real-time ultraviolet irradiance value at time t is given. Non-negativity constraints are applied to the calculated light-driven dynamic adjustment factor sequence, truncating any small negative values ​​caused by measurement noise to zero to ensure the thermodynamic rationality of the physical driving factors. The processed light-driven dynamic adjustment factor sequence is then structured and encapsulated according to timestamp indices to generate time-series weight modulation vectors corresponding to the directed edges of photosensitive oxidation types in the component degradation correlation graph. Through the above-mentioned quantification method of light-driven effects, the environmental stress monitoring data from the previous step is transformed into a light-driven dynamic adjustment factor sequence characterizing the intensity of photooxidation reaction kinetics. This achieves explicit guidance of the physical mechanism of the weights of photosensitive reaction paths in the graph neural network, improving the model's generalization ability and interpretability in predicting the stability of facial skincare products under varying light conditions.

[0078] S4.4: For directed edges in the component degradation correlation graph that have labels of acid-catalyzed hydrolysis, photosensitive oxidation and metal ion coordination acceleration, match the corresponding temperature dynamic regulation factor sequence, pH dynamic regulation factor sequence or light dynamic regulation factor sequence according to the edge type label, and perform multi-source factor weighted fusion processing to generate a comprehensive physical guidance edge weight matrix corresponding to each reaction path.

[0079] The algorithm iterates through all directed edges with physical semantic labels in the component degradation correlation graph, reads the type label attribute of each edge, and identifies its reaction mechanism category, including one or more combinations of acid-catalyzed hydrolysis, photosensitive oxidation, or metal ion coordination acceleration. For directed edges labeled as acid-catalyzed hydrolysis, values ​​matching the current time step and corresponding node environmental conditions are extracted from a pre-generated pH dynamic regulation factor sequence and used as the basic driving weights for this reaction path, reflecting the nonlinear regulatory effect of solution pH on the hydrolysis rate. For directed edges labeled as photosensitive oxidation, coefficient values ​​corresponding to real-time ultraviolet irradiance are obtained from the light dynamic regulation factor sequence, and combined with the thermal activation term in the temperature dynamic regulation factor sequence, a light-heat coupling driving weight is constructed to reflect the temperature dependence of photo-initiated free radical reactions.

[0080] For directed edges labeled as metal ion coordination acceleration types, the Arrhenius exponent term in the temperature dynamic adjustment factor sequence is extracted, and a pre-defined metal ion concentration correction coefficient is introduced to generate a comprehensive weight value reflecting the energy barrier reduction effect of coordination complex formation. For complex reaction edges simultaneously affected by multiple stress mechanisms, a weighted geometric average method is used to fuse the adjustment values ​​of each individual stress factor, and the comprehensive physical guidance edge weight is calculated using the following formula: in, This represents the combined physical guiding edge weight from node i to node j. , , These represent the dynamic adjustment factors of temperature, pH, and light, respectively, and k is the number of stress factors involved in the coupling.

[0081] The calculated integrated physical guide edge weights are filled into the corresponding positions in the adjacency matrix, replacing the original static topological connection weights, ensuring that the positions of non-zero elements in the matrix are strictly consistent with the directed edge structure of the component degradation correlation graph. Gradient truncation is performed on the generated integrated physical guide edge weight matrix to limit the upper and lower bounds of the weight values ​​to prevent numerical explosion or vanishing, ensuring numerical stability in the subsequent graph neural network message passing process. Through multi-source factor weighted fusion processing, the single-stress dynamic adjustment factor generated in the previous step is transformed into integrated physical guide edge weight matrix data characterizing the multi-stress coupling effect, achieving the expected technical effect of explicit quantitative mapping of graph topology to aging dynamics mechanism.

[0082] S4.5: Assign the generated integrated physical guidance edge weight matrix to the adjacency tensor structure of the component degradation correlation graph, replace the original static connection weights, and perform dynamic update operation of graph topology parameters to generate a parameterized differentiable physical guidance graph topology that supports backpropagation gradient flow and explicitly includes aging mechanism constraints.

[0083] Obtain the comprehensive physical guidance edge weight matrix generated in step S4.4. This matrix contains dynamic weight values ​​for each directed edge in the component degradation correlation graph, incorporating the multi-stress coupling effects of temperature, pH, and illumination, serving as the direct data source for updating the graph topology. Read the static adjacency tensor structure from the physical guidance graph topology structure output in step S3.5. This structure defines the connection relationships between nodes and the initialized fixed weight values, establishing the target data structure framework to be updated. Construct a differentiable weight placeholder tensor with dimensions completely identical to the static adjacency tensor structure. Initialize this placeholder using an automatic differentiation framework to ensure it can receive values ​​from the comprehensive physical guidance edge weight matrix and preserve gradient propagation paths.

[0084] A tensor masking operation is performed, using the Boolean mask of edge existence recorded in the physical guidance graph topology to precisely fill the non-zero dynamic weight values ​​in the comprehensive physical guidance edge weight matrix into the corresponding index positions of the differentiable weight placeholder tensor. For node pairs in the graph that are not connected by chemical reaction paths, the value is forcibly assigned to zero in the differentiable weight placeholder tensor to maintain the sparsity of the graph structure and avoid the consumption of ineffective computational resources and the introduction of noise. The filled dynamic weight tensor is then replaced element-wise with the original static adjacency tensor to generate a new parameterized adjacency tensor, in which each non-zero element is explicitly associated with the kinetic reaction rate coefficient under the current environmental stress conditions. A computational graph dependency relationship is established between the dynamic weight tensor and upstream symbolic regression parameters (such as activation energy, pre-exponential factor, and pH response coefficient) to ensure that when environmental stress parameters or symbolic rule parameters change slightly, the weight values ​​in the adjacency tensor can generate a corresponding gradient response through the chain rule. Perform topological consistency checks on the generated parameterized adjacency tensors to verify whether the weight values ​​of all active edges are within a physically reasonable range of positive real numbers, eliminate abnormal weights caused by numerical overflow or underflow, and ensure the numerical stability of the message passing process in the graph neural network.

[0085] Through the tensor mapping and dynamic replacement processing methods described above, the multi-stress coupled dynamic weights calculated in the previous step are transformed into a parameterized differentiable physical guide graph topology that supports backpropagation gradient flow, thereby realizing the explicit embedding and end-to-end optimization of aging mechanism constraints in the deep learning model.

[0086] Step S5: Construct the backbone of a graph neural network based on the parameterized differentiable physical guidance graph topology. Employ a hierarchical message passing mechanism to perform sparse multiplication of the adjacency matrix and gated aggregation operations. Use the linear combination coefficients generated in real-time by the symbolic rule parameters as message function constraints to output a sequence of latent vectors of component states explicitly guided by the aging mechanism. Specifically, this includes: S5.1: Sparsification is performed on the adjacency matrix in the parameterized differentiable physical guidance graph topology. Invalid connection edges are removed by using the reaction relationship labels between nodes to generate a highly sparsity component degradation association adjacency matrix, providing basic data support for reducing computational complexity in the future.

[0087] The system receives the parameterized differentiable physical guidance graph topology generated in step S4. This structure contains all active ingredient nodes, degradation product nodes, and a set of directed edges weighted according to kinetic rules. The adjacency tensor records the initial physical guidance edge weight matrix and the corresponding reaction type label. Each directed edge in the component degradation association graph is traversed, and its associated reaction type attribute set label, i.e., the reaction relationship label between nodes, is read. Specific mechanism categories such as acid-catalyzed hydrolysis, photosensitive oxidation, and metal ion coordination acceleration are identified, and the dynamic adjustment factor value of the edge at the current time step is extracted. A sparsity threshold judgment logic is set, and a non-zero validity check is performed on the comprehensive physical guidance edge weight of each directed edge. If the weight value is lower than the preset minimum reaction rate contribution threshold, the connection is determined to be an invalid chemical transformation path and marked as an object to be removed.

[0088] A high-sparseness mask matrix is ​​constructed, setting the corresponding positions of invalid edges in the adjacency matrix to zero, while retaining reaction path connections with significant physical meaning, thereby eliminating redundant topological connections caused by minor side reactions or noise interference. A sparsity reconstruction operation is performed on the adjacency matrix, storing the processed matrix in a compressed sparse row format, retaining only non-zero elements and their corresponding node indices and edge type labels, significantly reducing memory usage and computational load in subsequent message passing. Connectivity checks are performed on the reconstructed high-sparseness component degradation association adjacency matrix to ensure that each active component node has at least one valid input or output reaction path, eliminating completely isolated invalid nodes, and ensuring the physical integrity of the graph structure.

[0089] By using the above sparsification method, the dense physical guidance graph topology containing a large number of weak noise connections in the previous step is transformed into a highly sparse component degradation association adjacency matrix, which significantly reduces the computational complexity and explicitly highlights the key reaction paths, providing basic data support for the efficient message passing of the subsequent lightweight graph neural network.

[0090] S5.2: Based on the high sparsity component degradation association adjacency matrix, the symbolic rule parameters of each type of directed edge are extracted. Based on the dynamic adjustment factor mapping process, the temperature sensitive term, pH response term and light intensity coefficient are converted into real-time linear combination coefficients to generate a set of physical guidance weight coefficients for constraining the message passing process.

[0091] Based on the highly sparsity adjacency matrix of component degradation associations, a set of symbolic rule parameters corresponding to each type of directed edge is extracted. This set includes physical parameters such as reaction order, activation energy pre-exponential factor, and catalytic coefficient, which are resolved from the interpretable symbolic kinetic rule set generated in step S2. Each valid directed edge in the component degradation association graph is traversed, and its edge type label is identified, classifying the edge type into categories such as acid-catalyzed hydrolysis, photosensitive oxidation, metal ion coordination acceleration, or thermal degradation dominance. Dynamic regulation factor sequences matching each category are retrieved, including temperature dynamic regulation factors, pH dynamic regulation factors, and light dynamic regulation factors. For each directed edge, a linear combination coefficient mapping function is constructed. This function uses the symbolic rule parameters as the static basis and the dynamic regulation factors under real-time environmental stress as the modulation variables, calculating the linear combination coefficients at the current time step through weighted summation or product coupling. For reaction paths involving multi-stress coupling, a nonlinear activation function is used to fuse multi-source dynamic regulation factors, ensuring that the linear combination coefficients vary within the physically permissible range and avoiding gradient explosion or vanishing phenomena.

[0092] Specifically, for the first The node points to the first The linear combination coefficients of directed edges with nodes. The calculation follows the following mathematical logic: in, The basic bias term represents the intrinsic reaction rate constant under the condition of no external stress disturbance; The total number of environmental stress factors involved in the coupling; For the first The learnable weight coefficients of the stress-like factor are obtained by fine-tuning the parameters initialized by symbolic regression; For the first The dynamic adjustment factor function value corresponding to environmental stress. This is the environmental stress state vector corresponding to this edge.

[0093] If the edge is primarily affected by temperature, the environmental stress state vector corresponds to the current time value in the temperature dynamic adjustment factor sequence generated in step S4.1; if it is significantly affected by pH, it corresponds to the current time value in the pH dynamic adjustment factor sequence generated in step S4.2. All calculated linear combination coefficients are reorganized according to the graph's adjacency structure to generate a set of physical guidance weight coefficients with dimensions consistent with the high-sparseness adjacency matrix. Each element in this set explicitly carries physical mechanism constraint information under the current environmental conditions. Gradient truncation and normalization are performed on the physical guidance weight coefficient set to ensure that the coefficient distribution conforms to the non-negativity and order-of-magnitude characteristics of the chemical kinetic rate constant, preventing coefficient distortion caused by abnormal stress data.

[0094] Through the above dynamic mapping and coefficient generation processing, the static symbol rule parameters and real-time environmental stress dynamic adjustment factors obtained in the previous step are transformed into a set of physical guidance weight coefficients used to constrain the message passing process. This realizes the explicit following and dynamic adaptation of the graph neural network message passing mechanism to the aging physical mechanism, and significantly improves the prediction robustness and interpretability of the model under the varied accelerated aging experimental conditions.

[0095] S5.3: Perform hierarchical message passing operations on the initial embedding vector of the component state and the set of physical guiding weight coefficients. Use restricted linear combination as the core mechanism of the message function to dynamically fuse the state information of the upstream node with the edge weights to generate an intermediate layer message tensor containing local reaction dynamics features.

[0096] The initial embedding vector of the component state refers to a fixed-dimensional numerical vector generated for each node (representing an active ingredient or its degradation product) in the parameterized differentiable physics-guided graph topology, used to mathematically characterize the initial chemical state of the component. This initial embedding vector is directly derived from the node features in the graph structure constructed in the preceding step S4. Its specific value is initialized based on the physicochemical properties of the component (such as molecular weight, functional group type, initial concentration, etc.) and its corresponding symbol rule category (such as easy hydrolysis, photosensitivity, etc.), thereby encoding the component's domain knowledge into a mathematical form that can be processed by the graph neural network.

[0097] The initial embedding vector of the component state and the set of physical guidance weights are obtained as the input data basis for the hierarchical message passing operation. The initial embedding vector of the component state is initialized from the node features in the parameterized differentiable physical guidance graph topology generated in the preceding step S4, and the set of physical guidance weights is generated by mapping the symbolic rule parameters and dynamic adjustment factors of various types of directed edges extracted in S5.2. A message function of the form of a restricted linear combination is constructed, which is used to constrain the transmission process of upstream node state information to the current node. The core mechanism of the message function is defined as a linear transformation between the hidden state vector of the source node and the corresponding edge weights, ensuring that each step of information aggregation strictly follows the aging dynamics rules derived from symbolic regression. Each directed edge in the component degradation association graph is traversed, and the hidden state representation of the source node in the previous layer network is extracted. For any directed edge from node i to node j in the graph, the hidden layer feature vector of node i in the (l-1)th layer is obtained. This vector encodes the chemical stability state and local environmental response characteristics of the component at the previous time step.

[0098] The physical guidance weight coefficient associated with the directed edge is read. This coefficient is calculated in real time by dynamically adjusting the temperature sensitivity term, pH response term, and light intensity coefficient. The hidden state vector of the source node is multiplied element-wise or by matrix multiplication with the physical guidance weight coefficient to simulate the modulation effect of the chemical reaction rate on the amount of substance transformation.

[0099] The message vector passed from node i to node j is calculated using the following formula: in, This represents the message vector from node i to node j in the l-th layer message passing. This represents the physical guidance weight coefficient matrix determined by the edge type type(i,j). This matrix is ​​dynamically generated from the physical guidance weight coefficient set generated by S5.2. ⊙ represents the Hadamard product or linear projection operation. This represents the hidden layer feature vector of node i in the (l-1)th layer.

[0100] For multi-stress coupling scenarios, if the edge type involves the superposition of multiple reaction mechanisms, the weight coefficients of different types are weighted and fused. For example, for a degradation path simultaneously affected by heat and photo-oxidation, the heat-driven message component and the photo-driven message component are calculated separately, and then added together to obtain a comprehensive message vector to reflect the synergistic effect of multiple stresses. All incoming edge message vectors pointing to the target node j are initially aggregated to form an intermediate layer message tensor containing local reaction dynamics characteristics. This tensor retains the independent channel information of each upstream component's contribution to the degradation of the current component, without nonlinear activation or gating screening, to ensure the complete preservation of physical mechanism information.

[0101] By using a restricted linear combination message function processing method, the initial embedding vector of the component state and the set of physical guidance weight coefficients from the previous step are transformed into an intermediate layer message tensor containing local reaction dynamics characteristics. This enables explicit modeling and quantitative transmission of the chemical transformation relationship between components during the aging process of cosmetics, providing input data with clear physical meaning for subsequent gating aggregation operations.

[0102] S5.4: Perform gated aggregation operation based on the intermediate layer message tensor to generate aggregated node feature representations. Specifically, an adaptive gating unit is constructed using the sigmoid activation function to filter effective degradation path information from the intermediate layer message tensor, suppressing noise interference and retaining key response features.

[0103] The system receives an intermediate layer message tensor containing local reaction kinetics features generated in step S5.3. This tensor records the aggregated state information of neighboring nodes of each active component node at a specific time step, weighted by physically guided weighting coefficients. An adaptive gating unit based on the Sigmoid activation function is constructed. The intermediate layer message tensor is taken as input and mapped to a gating signal space with the same dimension as the node features through a linear transformation to quantify the confidence level of each degradation path's contribution to the current node's state update. The Sigmoid function is used to compress the linear output values ​​of the gating signal space to an open interval between 0 and 1, generating a gating coefficient matrix characterizing the strength of the physical consistency filtering. Values ​​close to 1 indicate that the path conforms to aging kinetics and the information is reliable, while values ​​close to 0 indicate that the path contains noise interference or violates physical constraints.

[0104] Where G is the generated gating coefficient matrix, σ is the Sigmoid activation function, and W g M is a learnable gated weight matrix. t Let b be the intermediate layer message tensor at time t. g This is a bias term.

[0105] The Hadamard Product operation is performed, and the generated gating coefficient matrix is ​​multiplied element-wise with the intermediate message tensor to achieve dynamic filtering of message content and suppress low-confidence information introduced by non-dominant reaction paths or experimental measurement errors.

[0106] Among them, M filtered This is the message tensor after physical consistency filtering, and ⊙ represents the Hadamard product operation.

[0107] The filtered message tensor is residually concatenated or spliced ​​with the previous hidden state of the current node to preserve the inherent chemical stability baseline information of the node and prevent the loss of intrinsic properties of key components due to over-filtering. Layer normalization is performed on the fused feature vector to eliminate numerical distribution shifts caused by differences in concentration magnitude between nodes of different components, ensuring gradient stability in subsequent nonlinear transformations. Through adaptive gating and Hadamard product filtering, the intermediate layer message tensor from the previous step is transformed into an aggregated node feature representation after physical consistency filtering. This achieves accurate preservation of effective degradation path information and efficient suppression of noise interference, significantly improving the model's robustness in representing complex multi-stress coupled aging behaviors.

[0108] S5.5: The aggregated node feature representation after physical consistency filtering is subjected to nonlinear transformation processing and mapped to the target latent space dimension through a fully connected layer to output the component state latent vector sequence explicitly guided by the aging mechanism, thus completing the forward inference process of the graph neural network backbone.

[0109] The system receives aggregated node feature representations filtered for physical consistency. These representations contain local response dynamics information filtered by gating units and serve as the input data basis for nonlinear transformation processing. A lightweight fully connected mapping layer is constructed, consisting of one-dimensional convolutional kernels or densely connected neurons. Its input dimension is consistent with the dimension of the aggregated node feature representation, and its output dimension is set to the target latent space dimension. This layer is used to achieve compressed mapping from high-dimensional aggregated features to low-dimensional component state latent vectors. Linear weighting is performed on the aggregated node feature representations. The weight matrix and bias vector of the fully connected layer are used to perform affine transformation on the input features to initially extract deep semantic features of components under specific aging environments, generating intermediate linear activation values. A nonlinear activation function is introduced to perform element-wise mapping processing on the intermediate linear activation values. The modified linear unit (ReLU) or exponential linear unit (ELU) is used to enhance the model's ability to fit complex degradation nonlinear relationships while avoiding the gradient vanishing problem, generating node latent features with nonlinear expressive capabilities.

[0110] The application layer normalization technique is used to standardize the latent features of nodes after nonlinear mapping, calculate the mean and variance of features within a batch, and adjust the feature distribution to zero mean and unit variance to accelerate model convergence and improve generalization adaptability to components of different concentration levels. The normalized features are then randomly sparsified using a discard regularization mechanism. During the training phase, some neuron outputs are randomly zeroed out with a preset probability, breaking the co-adaptive dependency between features and enhancing the robustness of the graph neural network backbone to noisy data. The regularized feature vectors are concatenated or stacked to form a state latent vector sequence containing all active components and their degradation product nodes. This sequence fully characterizes the internal state evolution of each component under the explicit guidance of the aging mechanism at the current time step.

[0111] Through nonlinear transformation and regularization of the fully connected layer, the aggregated node features obtained from the previous step after physical consistency filtering are transformed into a low-dimensional, dense, and highly generalizable latent vector sequence of component states. This achieves efficient encoding of component states during the complex aging process of cosmetics, providing high signal-to-noise ratio feature inputs for subsequent multi-step stability index prediction.

[0112] Step S6: Input the latent vector sequence of component states into a bi-objective loss function for backpropagation optimization to generate a graph neural network prediction model corrected for physical consistency. The bi-objective loss function includes: a primary task loss, used to calculate the multi-step prediction error of concentration, color, and viscosity stability indices; and an auxiliary task consistency regularization term, used to calculate the consistency regularization term of the sign rule parameters to force the evolution trajectory of edge weights to maintain consistency with the numerical trend of the dynamic expression. Specifically, it includes: S6.1: Perform fully connected mapping on the latent vector sequence of component states to generate a multi-step prediction result set of stability indices containing concentration prediction, colorimetric prediction, and viscosity prediction, which serves as the direct input object for the main task error calculation.

[0113] S6.2: Perform mean squared error calculation based on the multi-step prediction result set of stability index and real accelerated aging experimental data to generate a scalar of the main task prediction error that characterizes the model fitting accuracy, which is used to guide the gradient update direction of the model's basic prediction capability.

[0114] S6.3: Utilize the temperature-sensitive term coefficients in the interpretable symbolic dynamics rule set and the directional cosine similarity calculation of the directed edge weight evolution trajectory in the current training round to generate symbolic rule parameter consistency regularization term values ​​that characterize the degree of adherence to physical mechanisms.

[0115] The real-time evolution trajectory vectors of the directed edge weights of the lightweight graph neural network in the current training epoch are extracted as a sequence of physical state observations representing the dynamics of data-driven learning. This trajectory vector records the numerical changes in the edge weights representing specific chemical reaction paths as they are updated with gradient descent from the initial time step to the current time step, reflecting the model's dynamic fitting process to the degradation rate. The temperature-sensitive term coefficient vectors corresponding to the reaction type in the interpretable symbolic kinetic rule set are analyzed. These coefficient vectors are obtained by linearizing the activation energy and pre-exponential factor in the Arrhenius equation, representing the theoretical sensitivity benchmark of the reaction rate to temperature changes under ideal physical constraints. Standardization preprocessing is performed on the edge weight evolution trajectory vectors and the temperature-sensitive term coefficient vectors to eliminate dimensional differences and map them to the same feature space, ensuring that the geometric meaning of subsequent similarity calculations is clear and unaffected by numerical scale.

[0116] A cosine similarity calculation model is constructed, treating the standardized edge weight evolution trajectory as a direction vector in a high-dimensional space, and the temperature-sensitive term coefficient as a reference direction vector. The cosine of the angle between the two vectors is calculated to quantify the consistency between data-driven trends and physical prior trends. The following formula is used to calculate the... The reaction pathway is in the first Directional cosine similarity across training rounds: in, For the first The weight evolution trajectory vector of the edge in the current training round. For the first The temperature-sensitive coefficient vector corresponding to each edge. This represents the vector dot product operation. The L2 norm of a vector. For the first The reaction pathway is in the first Directional cosine similarity across training rounds.

[0117] A global average aggregation operation is performed on the direction cosine similarity values ​​of all valid reaction paths in the graph to generate a scalar index characterizing the overall physical compliance of the model. This index reflects the directional consistency between the degradation kinetics learned by the neural network and the principles of chemical thermodynamics. A non-negative constraint mechanism is introduced to map the direction cosine similarity to the interval [0, 1]. When the similarity is close to 1, it indicates that the data-driven trend is highly consistent with the physical mechanism, while when the similarity is close to 0 or negative, it indicates a serious physical violation. Based on the mapped similarity index, a consistency regularization loss function is constructed. By inverting or reciprocalizing the function, the higher the physical consistency, the smaller the regularization loss, thereby guiding the model parameters to optimize in a direction consistent with physical laws during backpropagation.

[0118] S6.4: Perform a weighted summation operation on the scalar of the main task prediction error and the consistency regularization term of the sign rule parameter to generate a dual-objective joint loss function value that integrates data-driven accuracy and physical interpretability constraints, which serves as the unified objective function for backpropagation optimization.

[0119] The system receives the scalar value of the main task prediction error and the value of the consistency regularization term of the sign rule parameter generated in the preceding steps, serving as the input data source for constructing the bi-objective joint loss function. The scalar value of the main task prediction error characterizes the model's fitting deviation on macroscopic stability indicators such as concentration, color, and viscosity, reflecting the prediction accuracy at the data-driven level. The value of the consistency regularization term of the sign rule parameter characterizes the degree of deviation between the evolution trajectory of the graph neural network edge weights and the trend of the physical dynamics expression, reflecting the constraint compliance at the physical mechanism level. A hyperparameter balancing coefficient is set to adjust the relative weights of data-driven accuracy and physical interpretability constraints in the total loss function. This balancing coefficient is dynamically configured based on the signal-to-noise ratio of the accelerated aging experimental data and the confidence level of the physical prior knowledge. Typically, physical constraints are given higher weights in the early stages of training to guide the model to converge to the physical feasible region, and the weights of physical constraints are gradually reduced in the later stages of training to fine-tune the prediction accuracy.

[0120] A weighted summation operation is performed, multiplying the scalar of the main task's prediction error by the data-driven weight coefficient, and multiplying the value of the consistency regularization term of the sign rule parameter by the physical constraint weight coefficient. The two loss components are then fused using a linear superposition method to eliminate the impact of dimensional differences on gradient updates, ensuring that both tasks have a balanced gradient contribution rate during backpropagation. The joint loss function value for the two objectives is calculated using the following formula: Where L is the value of the joint loss function for the two objectives; λ1 is the data-driven weighting coefficient, used to control the degree of influence of the prediction error of the main task; L MSE λ is the scalar of the prediction error for the main task, i.e., the mean squared error loss; λ² is the physical constraint weighting coefficient, used to control the influence of the consistency regularization term of the sign rule parameters; L Phys This is the value of the consistency regularization term for the symbol rule parameters, which is the complement of the direction cosine similarity loss.

[0121] The calculated joint loss function value of the two objectives is scalarized to generate a backpropagation target signal with a single gradient. This signal simultaneously contains the gradient direction that improves prediction accuracy and the gradient direction that enhances physical consistency, serving as a unified optimization objective for subsequent adaptive moment estimation algorithms.

[0122] S6.5: An adaptive moment estimation algorithm is executed based on the joint loss function value of the two objectives to synchronously update the linear combination coefficients of the message function of the graph neural network backbone and the parameters of the dynamic adjustment factor, and finally generate a graph neural network prediction model after physical consistency correction.

[0123] The system receives the joint loss function value for both objectives from step S6.4. This value integrates the scalar of the main task prediction error and the value of the consistency regularization term for the sign rule parameters, serving as the global optimization objective for the current training epoch. A gradient update mechanism based on Adaptive Moment Estimation (Adam) is constructed, initializing the first-order and second-order moment estimation variables as zero vectors. A learning rate decay strategy and momentum hyperparameter are set to provide an initial state baseline for subsequent parameter iterations.

[0124] Gradient calculation is performed on the linear combination coefficients of the message function in the backbone of the graph neural network. The bi-objective joint loss function value is backpropagated using the chain rule. The partial derivatives of the loss function with respect to the temperature-sensitive term coefficient, pH response term coefficient, and light intensity coefficient are calculated to generate the original gradient tensors corresponding to each physical driving factor. Synchronous gradient calculation is performed on the dynamic adjustment factor parameters, tracing the complete computational graph path from the edge weight dynamic adjustment factor to the adjacency matrix sparse multiplication, then to the gated aggregation operation and the final stability index prediction value. This obtains the sensitivity information of the dynamic adjustment factor relative to the joint loss function, forming a parameter update vector containing the physical constraint gradient. The first-order moment estimate is updated using the exponentially weighted moving average method. The original gradient tensor calculated at the current time step is weighted and fused with the first-order moment estimate from the previous time step to eliminate gradient oscillation noise, extract the trend direction of gradient change, and generate a smoothed first-order gradient estimate.

[0125] The second-order moment estimate is updated using an exponentially weighted moving average method with squared gradients. The element-wise squared value of the original gradient tensor at the current time step is weighted and fused with the second-order moment estimate from the previous time step to capture the fluctuation range of the gradient magnitude, generating a second-order gradient estimate for adaptive learning rate adjustment. Bias correction is performed on the first-order and second-order moment estimates. Correction coefficients are calculated based on the current training steps to eliminate the bias towards zero in the initialization phase of the moment estimates, generating unbiased corrected first-order and second-order gradient estimates, ensuring update stability in the early training phase. An adaptive learning rate is calculated based on the bias-corrected moment estimates. The unbiased corrected first-order gradient estimate is divided by the square root of the unbiased corrected second-order gradient estimate, with a small constant added to prevent division by zero errors. This enables differentiated step size adjustment for different parameter dimensions, generating a dedicated adaptive update step size for each parameter to be updated. The generated adaptive update step size is applied to the linear combination coefficients of the message function and the dynamic adjustment factor parameters. The parameter subtraction update operation is performed, and the network weights are adjusted along the negative gradient direction. This allows the model to minimize the prediction error while forcing the evolution trajectory of the edge weights to conform to the trend constraint of the symbolic dynamics rule.

[0126] Step S7: For components not seen in the new formulation, reuse the embedding space of existing component nodes, and fine-tune the sign rule adaptation coefficients between them and adjacent nodes using a small number of samples to generate an adaptive prediction model instance with cross-formula transfer capabilities. Specifically, this includes: S7.1: Obtain the list of unseen ingredients in the new formula and the corresponding small amount of accelerated aging experimental concentration time-series data. Perform semantic similarity matching operation based on the standardized graph node set generated in the previous step, and map the unseen ingredients to the existing component node embedding space with the closest chemical structure or functional properties to generate an unseen component mapping result set containing the initial latent vector representation.

[0127] S7.2: Based on the unseen component mapping result set and the interpretable symbolic dynamics rule set output by the previous steps, extract the potential reaction path information between known component nodes adjacent to the unseen component, and use linear interpolation to initialize the initial value of the symbolic rule adaptation coefficient of the edge weight between the unseen component and the adjacent node to generate the initial matrix of symbolic rule adaptation coefficients to be optimized.

[0128] Based on the unseen component mapping result set generated by S7.1, the K known component neighbor nodes with the highest semantic similarity in the normalized graph node set for each unseen component node are extracted to construct a local neighborhood topological subgraph. This step aims to utilize the physicochemical properties of known components as a priori benchmark to establish an initial reaction kinetic context for unseen components.

[0129] The interpretable symbolic kinetic rules set output from step S2 is retrieved, and all reaction path rules involving the aforementioned K known neighbor nodes are selected. For each selected reaction path, its corresponding symbolic kinetic expression structure is analyzed to identify key physical parameters determining the reaction rate, including the Arrhenius activation energy, pre-exponential factor, pH dependence index, and photoluminescence quantum yield coefficient. The Euclidean distance between the molecular descriptor vectors of the unseen component nodes and the molecular descriptor vectors of the K known neighbor nodes is calculated to generate a normalized structural similarity weight vector. This weight vector characterizes the proximity of the unseen component to each known neighbor in the chemical structure space, serving as the basis for confidence in subsequent linear interpolation. For potential connection edges between the unseen component and the i-th known neighbor node, the optimized symbolic rule fitting coefficient vector of that neighbor node in historical training data is extracted. Using linear interpolation with the structural similarity weights as interpolation coefficients, a weighted average is performed on the fitting coefficient vectors of the K neighbor nodes to generate an initial estimate of the symbolic rule fitting coefficients for the unseen component.

[0130] Specifically, for the unseen component u and the known component v i The edge weight adaptation coefficients between them are initialized using the following linear interpolation formula: Among them, w j This represents the relationship between the unseen component u and its j-th known neighbor node v. j Normalized structural similarity weights between them; Indicates a known neighbor node v j With the target known component v i The sign rule adaptation coefficients determined in the pre-trained model; K is the number of nearest neighbor nodes selected.

[0131] If multiple reaction type labels exist between an unseen component and a known neighbor node (e.g., simultaneous acid-catalyzed hydrolysis and photosensitive oxidation), the linear interpolation process described above is performed independently for each reaction type to generate a multi-dimensional initial value vector of fitting coefficients. This ensures differentiated initialization of the kinetic response characteristics under different stress conditions. All calculated initial sign rule fitting coefficients are assembled into a sparse matrix, where the row indices correspond to unseen component nodes, and the column indices correspond to their adjacent known component nodes and reaction type combinations. Non-negativity constraints are applied to this matrix, truncating any negative coefficients caused by interpolation errors to zero to ensure the physical rationality of the kinetic rate constant.

[0132] S7.3: Perform small-sample gradient descent optimization on the initial matrix of the symbolic rule fitting coefficients. Use a small amount of experimental concentration time-series data of the new formulation as a supervision signal to calculate the prediction error of the main task and backpropagate to update the fitting coefficients. At the same time, freeze the embedding parameters and core message passing mechanism of the existing component nodes to generate the optimized matrix of the symbolic rule fitting coefficients after fine-tuning.

[0133] The initial matrix of adaptation coefficients for the symbol rules to be optimized, generated in step S7.2, is received and loaded into the computational graph memory as a trainable parameter tensor. Simultaneously, the gradient update flags for the embedding vectors of all nodes, the weights of the message passing layer, and the parameters of the gating units in the pre-trained graph neural network backbone are locked to ensure that parameter iteration is performed only for the edge weight adaptation coefficients introduced by the new formulation during subsequent backpropagation. A few-sample supervised learning loop based on a small amount of accelerated aging experimental data of the new formulation is constructed. The initial latent vector representation of the unseen components obtained in step S7.1 is concatenated with the node features of the known components and input into the frozen graph neural network backbone. Forward inference is performed to obtain the predicted sequence of stability indicators under the current adaptation coefficients.

[0134] Define the main task loss function for the small sample fine-tuning stage, use the mean square error to measure the deviation between the predicted concentration time series curve and the actual HPLC detection data, and calculate the main task loss value using the following formula: Where N is the number of component nodes participating in fine-tuning, and T is the total number of time steps. Let be the actual concentration observation value of the i-th component at time t. To predict concentration values ​​for the model.

[0135] Introducing a physical consistency regularization constraint, we calculate the cosine similarity penalty between the direction of the derivative of the sign rule adaptation coefficient in the current iteration and the evolution direction of the standard dynamic rule parameters determined in step S6. An auxiliary loss term is constructed using the following formula: Among them, wadapt To adapt the coefficient vector for the current sign rule to be optimized, w std The standard kinetic rule parameter vector is obtained from the migration of known components. This constraint term forces the degradation rate of the new components to conform to the basic chemical kinetic laws.

[0136] The main task loss and regularization loss are combined to form a joint optimization objective function. Regularization weights are set to balance data fitting accuracy and physical interpretability. An automatic differentiation engine is used to calculate the gradient of the joint loss function relative to the initial matrix of the symbolic rule fitting coefficients. Since the backbone network parameters are frozen, the gradient flow only passes through local computational paths related to the new formula's edge weights. Mini-batch gradient descent is executed, updating the symbolic rule fitting coefficient matrix based on the calculated gradient values. An adaptive learning rate scheduling strategy is adopted, using a larger step size in the early stages of training to quickly approach the optimal solution region, and gradually decreasing the step size in later stages to finely adjust the coefficient values ​​and avoid overfitting on small sample data. The prediction error on the validation set is monitored. When the main task loss no longer decreases significantly or the validation error rises within several consecutive iterations, an early stopping mechanism is triggered, terminating the gradient update process and saving the current optimal state of the symbolic rule fitting coefficient matrix.

[0137] S7.4: Based on the symbol rule adaptation coefficient optimization matrix, update the directed edge weight dynamic adjustment factor involving unseen components in the component degradation correlation graph, substitute the optimized adaptation coefficient into the parameterized differentiable physical guidance graph topology constructed in the previous step, and perform a local topology reconstruction operation to generate a local physical guidance graph topology instance that adapts to the characteristics of the new formulation.

[0138] The system receives a finely tuned and corrected symbolic rule adaptation coefficient optimization matrix. This matrix contains dynamic adjustment parameters for the reaction paths between unseen components and adjacent known component nodes in the new formulation, serving as the core input data for local topology reconstruction. It traverses all directed edges involving unseen components in the component degradation correlation graph, extracting the corresponding reaction type label for each edge, including physical semantic categories such as acid-catalyzed hydrolysis, photosensitive oxidation, or metal ion coordination acceleration, to determine the computational logic branch for subsequent dynamic adjustment factors. Based on the reaction type label, it retrieves matching kinetic expression prototypes from the interpretable symbolic kinetic rule set, obtains the positions of the symbolic rule adaptation coefficients to be updated in the expression, and fills the corresponding positions with the values ​​from the symbolic rule adaptation coefficient optimization matrix, replacing the original general coefficients or initial interpolation coefficients. It reads the real-time temperature, pH, and light intensity values ​​from the current accelerated aging experimental environmental stress parameter labels, substitutes these environmental stress parameters into the updated kinetic expression prototype, performs multi-stress coupling effect calculations, and generates a sequence of edge weight dynamic adjustment factors for a specific time step.

[0139] For the temperature-sensitive term, the calculation is performed using the Arrhenius equation, and the temperature dynamic adjustment factor is determined using the following formula: in, As a temperature dynamic adjustment factor, and These are the optimized pre-exponential factor and activation energy fit coefficient, respectively. Let be the ideal gas constant. The value is the real-time absolute temperature. For the pH response term, the pH dynamic adjustment factor is calculated using the following formula, based on the acid-catalyzed or base-catalyzed rate equation model: in, As a pH dynamic regulator, The optimized acid catalytic rate constant is the fitting coefficient. This represents the real-time pH value. For the light intensity coupling term, based on the relationship between photon yield and the product of light intensity, the light dynamic adjustment factor is calculated using the following formula: in, As a dynamic adjustment factor for illumination, The optimized quantum yield adaptation coefficient, This represents real-time ultraviolet irradiance.

[0140] The calculated dynamic adjustment factors for temperature, pH, and illumination are weighted and fused, and weights are assigned according to the primary and secondary relationships of the reaction mechanisms to generate a comprehensive physical guidance edge weight value for each directed edge involving unseen components. This comprehensive physical guidance edge weight value is assigned to the corresponding position in the adjacency tensor structure of the component degradation correlation graph, replacing the original static connection weights or the weight values ​​from the previous iteration, thus achieving local dynamic updates of the graph topology parameters. A sparsity check is performed on the updated adjacency tensor structure to ensure that the original edge weights not involving unseen components remain unchanged, and only the topology of the local subgraphs related to the new components is reconstructed, maintaining the stability of the overall graph structure and computational efficiency.

[0141] S7.5: Integrate the local physical guidance graph topology instance with the unchanged original graph neural network backbone parameters to encapsulate a complete adaptive prediction model instance. This instance retains the physical interpretability constraints of the original model and has the ability to respond quickly to new ingredients, so as to output an adaptive prediction model instance with cross-formula transfer capability for subsequent aging process simulation.

[0142] Step S8: Utilize the adaptive prediction model instance to extrapolate the accelerated aging process of the current cosmetic product, outputting predicted stability indices and component-level attribution analysis results for future time steps, thus completing the lightweight graph neural network prediction process guided by aging dynamics symbolic regression. Specifically, this includes: S8.1: Obtain the initial component composition vector of the current cosmetic formula to be tested and the preset environmental stress parameter sequence of accelerated aging experiment. Based on the cross-formula migration strategy in the adaptive prediction model instance, map the unseen components to the embedding space of the existing component nodes and fine-tune the symbol rule adaptation coefficient to generate an initial physical guidance graph topology structure with specific formula adaptability.

[0143] S8.2: Based on the initial physical guidance graph topology structure, a hierarchical message passing mechanism is executed. The linear combination coefficients generated in real time by the symbol rule parameters are used to constrain the sparse multiplication and gated aggregation operations of the adjacency matrix. The state hidden vectors of each active ingredient and its degradation product node in the graph are iteratively updated to generate a sequence of component state hidden vectors that reflect the multi-step aging evolution trajectory.

[0144] The system receives an initial physical guidance graph topology with specific recipe adaptability and extracts its set of initial state latent vectors for nodes, a directed edge adjacency matrix with physical semantic labels, and a sequence of dynamically adjusted factors generated in real time by symbolic rule parameters. These serve as the input data basis for a hierarchical message passing mechanism. The adjacency matrix in the initial physical guidance graph topology is subjected to sparsification encoding. Invalid and zero-weight edges are removed based on the reaction relationship labels between nodes, constructing a highly sparse component degradation-related adjacency matrix, significantly reducing the computational complexity and memory consumption of subsequent matrix multiplication operations.

[0145] Based on a highly sparsity adjacency matrix for component degradation association, the algorithm traverses each directed edge in the graph labeled with acid-catalyzed hydrolysis, photosensitive oxidation, or metal ion coordination acceleration. It then invokes a dynamic adjustment factor mapping process, substituting the environmental stress parameters of the current time step into the corresponding symbolic kinetic rule expression to calculate real-time linear combination coefficients. Using restricted linear combination as the core mechanism of the message function, it performs element-wise weighted operations on the state latent vector of the upstream source node and the real-time linear combination coefficients of the corresponding edges, fusing local reaction kinetic features to generate an intermediate layer message tensor characterizing the potential for single-step chemical transformation. For the multiple incoming edge message tensors converged at the target node, a gated aggregation operation is performed. An adaptive gating unit is constructed using the Sigmoid activation function, filtering effective degradation path information based on the physical confidence of each reaction path, suppressing noise interference, and retaining key reaction features to generate a physically consistent, aggregated node feature representation.

[0146] The aggregated node feature representations, after physical consistency filtering, undergo nonlinear transformation processing. These representations are then mapped to the target latent space dimension via a fully connected layer. This updates the state latent vectors of each active ingredient and its degradation product node in the graph, completing a single-layer message passing iteration. This hierarchical message passing, gated aggregation, and state update process is repeated until a preset time step or convergence threshold is reached. The final output is a sequence of component state latent vectors reflecting the multi-step aging evolution trajectory.

[0147] S8.3: Perform decoding mapping processing on the latent vector sequence of the component state, convert the high-dimensional latent space features into physically interpretable stability index values ​​through the model output layer, calculate the residual concentration of active ingredients, the change value of system color and the fluctuation of rheological viscosity at multiple future time steps, so as to generate a set of multi-step stability index prediction values.

[0148] The component state latent vector sequence output from step S8.2 is received. This sequence contains high-dimensional feature representations of each active component node and its degradation product node in the graph at the current time step and the predicted future time step, serving as the input data source for decoding and mapping. A multi-task parallel decoder architecture is constructed. For three types of physically interpretable stability indicators—concentration residual rate, system color change value, and rheological viscosity fluctuation—independent fully connected neural network branches are instantiated to ensure decoupling mapping of features of different physical dimensions. A global pooling operation is performed on the component state latent vector sequence, and a weighted averaging mechanism is used to aggregate the feature information of all nodes in the graph to generate a context-aware vector representing the overall formulation state, eliminating the influence of node arrangement order on the prediction results.

[0149] The context-aware vector is input to the concentration prediction branch. Through a combination of two linear transformations and a nonlinear activation function, it is mapped to an output dimension that matches the number of active ingredients, directly representing the molar concentration or mass fraction of each key component at a specific time point. A physical constraint layer is introduced to correct the output of the concentration prediction branch. A non-negative activation function is used to force all predicted concentration values ​​to be greater than zero, and a sum-conservation constraint is applied to ensure that the sum of the concentrations of the parent component and the main degradation products does not exceed the initial total feed amount. The residual concentration rate of active ingredients is calculated using the following formula: in, Let be the residual concentration of the i-th active ingredient. Let be the concentration of the i-th component predicted by the model at time t. Let be the initial concentration of the i-th component.

[0150] The context-aware vector is input into the chromaticity prediction branch, and combined with the chromaticity reference value at historical moments, the Lab color space coordinate offset of the future time step is output through regression mapping, focusing on capturing the yellowing trend caused by oxidation or Maillard reaction.

[0151] The colorimetric change value of the system is calculated using the colorimetric difference formula to quantify the degree of degradation of the formulation's appearance stability. in, This is the total color difference value. , , These represent the coordinate differences between the predicted time and the initial time on the brightness, red-green axis, and yellow-blue axis, respectively.

[0152] The context-aware vector is input into the viscosity prediction branch to extract the implicit rheological structure factor. The dynamic viscosity coefficient is output through nonlinear fitting, which reflects the texture change caused by the destruction of the emulsion system or the degradation of polymers.

[0153] Calculate the rheological viscosity fluctuation to assess the degree of deviation from the physical stability of the system: in, This represents the relative fluctuation in viscosity. To predict the dynamic viscosity at a given time, This represents the initial dynamic viscosity.

[0154] The original predicted values ​​output from the above three branches are subjected to inverse normalization. Based on the statistical distribution parameters of each indicator recorded during the training phase, the standardized values ​​are restored to engineering unit values ​​with actual physical meaning.

[0155] The concentration residual rate, color change value and viscosity fluctuation amount of multiple time steps are integrated and reorganized into a structured data matrix according to the time index order to form a set of multi-step stability index prediction values ​​covering the entire prediction period.

[0156] S8.4: Based on the state change gradient and edge weight evolution trajectory of each node in the component state latent vector sequence, backtrack analysis is performed on the contribution of different environmental stress factors to specific chemical reaction paths, and key degradation nodes and dominant reaction types that lead to the deterioration of stability indicators are identified, so as to generate component-level attribution analysis results that include the location of key failure components and mechanism attribution.

[0157] The component state latent vector sequence output from step S8.2 and the dynamic edge weight evolution trajectory in the topology of the parameterized differentiable physical guide graph generated in step S4 are obtained as the input data basis for attribution analysis.

[0158] Gradient calculations are performed on the node features at each time step in the latent vector sequence of component states. Automatic differentiation techniques are used to solve for the partial derivatives of the target stability index with respect to the latent states of each active component node, quantifying the instantaneous sensitivity of each component's concentration changes to overall stability. Based on the calculated node state gradients, combined with the graph adjacency matrix and edge weights of the current time step, a reverse message passing operation is performed to propagate the error signal from the output layer back along directed edges to the input layer nodes, tracing the upstream critical reaction paths leading to stability degradation. The time-evolving sequence of dynamic adjustment factors recorded in step S4, including temperature sensitivity, pH response, and light intensity coefficients, is extracted. These are then multiplied element-wise with the path contribution obtained from the reverse propagation to separate the contribution weights of different environmental stress factors to specific chemical reaction rates.

[0159] A stress-path coupling contribution matrix is ​​constructed, mapping the reaction type label corresponding to each directed edge to its associated environmental stress factor weight. By accumulating the contributions of all edges under the same reaction type, a global influence score for various degradation mechanisms (such as acid-catalyzed hydrolysis and photosensitive oxidation) is generated. Reaction paths with global influence scores exceeding a preset threshold and their starting nodes are identified. These nodes are marked as key failure components, and the corresponding high-weight reaction types are marked as dominant degradation mechanisms, forming a preliminary attribution candidate set. A temporal consistency check is performed on the key failure components in the attribution candidate set, checking whether their state gradient signs continuously point to the direction of concentration decrease within consecutive time steps. Nodes with transiently abnormally brightened states due to numerical noise are removed to ensure the physical rationality of the attribution results. The validated list of key failure components, dominant degradation mechanism types, and specific contribution values ​​of each stress factor are integrated to generate structured component-level attribution analysis results, clearly indicating the core chemical mechanisms and environmental triggers leading to decreased cosmetic stability.

[0160] S8.5: Integrate the predicted values ​​of the multi-step stability index with the results of the component-level attribution analysis to construct a visual evaluation report that includes a time-dimensional trend curve and a component contribution heatmap, and output the final aging process deduction conclusion to complete the entire process of lightweight graph neural network prediction based on aging dynamics symbolic regression.

[0161] For those skilled in the art, various other corresponding changes and modifications can be made based on the technical solutions and concepts described above, and all such changes and modifications should fall within the protection scope of the claims of this invention.

[0162] Unless otherwise defined, the technical or scientific terms used herein shall have the ordinary meaning as understood by one of ordinary skill in the art to which this application pertains. The terms “first,” “second,” “third,” and similar terms used in this patent application specification and claims do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Similarly, the terms “an” or “a” and similar terms do not indicate a quantity limitation, but rather indicate the presence of at least one. The terms “comprising” or “including” and similar terms mean that the elements or objects preceding “comprising” or “including” encompass the elements or objects listed following “comprising” or “including” and their equivalents, and do not exclude other elements or objects. The “multiple” mentioned in the embodiments of this application refers to two or more. A and / or B indicate three possibilities: A; B; and A and B.

[0163] The above description is merely an exemplary embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and such modifications or substitutions should all be covered within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for predicting the time-series trend of accelerated aging data for cosmetic stability testing, characterized in that, include: S1: Obtain the time-series concentration data of each active ingredient in the typical cosmetic ingredient library under multiple stress conditions in accelerated aging experiments, and record the corresponding environmental stress parameter labels to form the original multidimensional dataset; S2: Perform symbolic regression on the original multidimensional dataset to generate an interpretable symbolic dynamics rule set; S3: Using the interpretable symbolic dynamics rule set to define the graph topology, generate a component degradation correlation graph; S4: Calculate the dynamic adjustment factor according to the interpretable symbolic dynamic rule set, and assign the dynamic adjustment factor to the directed edge weight of the corresponding type label in the component degradation association graph to generate the physical guidance graph topology. S5: Construct the backbone of the graph neural network based on the physical guidance graph topology, use a hierarchical message passing mechanism to perform sparse multiplication of the adjacency matrix and gated aggregation operations, and convert the symbolic rule parameters into linear combination coefficients to output the component state hidden vector sequence; S6: Input the latent vector sequence of the component states into the bi-objective loss function for backpropagation optimization to generate a graph neural network prediction model; S7: For components not seen in the new formula, reuse the embedding space of existing component nodes, fine-tune the sign rule adaptation coefficients between them and adjacent nodes through a small number of samples, and generate adaptive prediction model instances. S8: Using the aforementioned adaptive prediction model instance, the accelerated aging process of the current cosmetics is simulated, and the predicted values ​​of stability indicators and component-level attribution analysis results for future time steps are output.

2. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 1, characterized in that, The graph topology is defined using the interpretable symbolic dynamics rule set, including: mapping each active ingredient and its key degradation products in the typical cosmetic ingredient library as graph nodes, and constructing directed edges with reaction type labels based on the reaction relationships identified in the interpretable symbolic dynamics rule set.

3. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 2, characterized in that, The reaction types include acid-catalyzed hydrolysis, photosensitive oxidation, and metal ion coordination acceleration.

4. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 3, characterized in that, Step S3 specifically includes: Based on the list of active ingredients and their key degradation products identified in the set of interpretable symbolic kinetic rules, a unique identifier mapping process is performed on each chemical component entity to generate a standardized graph node set containing specific substance instances. By utilizing the chemical transformation logic between nodes in the standardized graph node set, and based on the reaction path described by the interpretable symbolic kinetic rule set, a directed connection relationship determination operation is performed to generate a potential edge connection matrix. For each valid connection edge in the potential edge connection matrix, the reaction mechanism description information corresponding to the interpretable symbolic kinetic rule set is extracted, and the classification and labeling processing of acid-catalyzed hydrolysis, photosensitive oxidation or metal ion coordination acceleration is performed to generate a reaction type attribute set. By combining the standardized graph node set, potential edge connection matrix, and reaction type attribute set, graph structure assembly and topology integrity verification operations are performed to generate a component degradation correlation graph. Based on the generated component degradation correlation graph, all directed edges are traversed, the category information in the reaction type attribute set is solidified into edge metadata, and the physical guidance graph topology is output.

5. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 4, characterized in that, Specific examples of the substances mentioned include: nicotinamide and vitamin E acetate.

6. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 4, characterized in that, Step S4 specifically includes: Based on the temperature-sensitive terms parsed from the set of interpretable symbolic dynamic rules, the Arrhenius activation energy parameter and pre-exponential factor are extracted. Combined with the real-time temperature values ​​in the environmental stress parameter labels of the current accelerated aging experiment concentration time series data, an exponential function operation is performed to generate a temperature dynamic adjustment factor sequence. Using the pH response terms identified in the set of interpretable symbolic kinetic rules, the real-time pH values ​​in the environmental stress parameter labels are read, and then substituted into the preset acid catalysis or base catalysis rate equation model for nonlinear mapping calculation to generate a pH dynamic adjustment factor sequence. Based on the light intensity coefficient defined in the interpretable symbolic dynamics rule set, real-time ultraviolet irradiance data in the environmental stress parameter label is obtained, and the photon yield and light intensity product operation is performed to generate a light dynamic adjustment factor sequence. For the directed edges in the component degradation correlation graph that have labels of acid-catalyzed hydrolysis, photosensitive oxidation and metal ion coordination acceleration, the corresponding temperature dynamic regulation factor sequence, pH dynamic regulation factor sequence or light dynamic regulation factor sequence are matched according to their edge type labels, and multi-source factor weighted fusion processing is performed to generate a comprehensive physical guided edge weight matrix. The generated integrated physical guidance edge weight matrix is ​​assigned to the adjacency tensor structure of the component degradation correlation graph, replacing the original static connection weights. Dynamic update operation of graph topology parameters is then performed to generate the physical guidance graph topology structure.

7. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 6, characterized in that, Step S5 specifically includes: The adjacency matrix in the physical guidance graph topology is sparsified, and invalid connection edges are removed using the reaction relationship labels between nodes to generate a component degradation association adjacency matrix. Based on the component degradation association adjacency matrix, the symbolic rule parameters of each type of directed edge are extracted. Based on the dynamic adjustment factor mapping process, the temperature sensitive term, pH response term and light intensity coefficient are converted into real-time linear combination coefficients to generate a set of physical guided weight coefficients. A hierarchical message passing operation is performed on the initial embedding vector of the component state and the set of physical guidance weight coefficients. A restricted linear combination is used as the core mechanism of the message function to dynamically fuse the state information of the upstream node with the edge weights to generate the intermediate layer message tensor. Gated aggregation operations are performed based on intermediate layer message tensors to generate aggregated node feature representations. The aggregated node feature representations after physical consistency filtering are subjected to nonlinear transformation processing and mapped to the target latent space dimension through a fully connected layer, outputting the component state latent vector sequence.

8. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 7, characterized in that, The gated aggregation operation based on the intermediate layer message tensor includes: constructing an adaptive gating unit using the sigmoid activation function, filtering effective degradation path information on the intermediate layer message tensor, suppressing noise interference, and retaining key reaction features.

9. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 1, characterized in that, The dual-objective loss function includes: a primary task loss, used to calculate the multi-step prediction error of concentration, color, and viscosity stability indices; and an auxiliary task consistency regularization term, used to calculate the consistency regularization term of the sign rule parameters.

10. The method for predicting the time-series trend of accelerated aging data for cosmetic stability testing according to claim 1, characterized in that, The multi-stress conditions include: heat, light, oxygen, and pH.

Images