Regulation and control method and system for environment and equipment cluster, electronic equipment and storage medium

Abstract

The invention provides an environment and equipment cluster regulation and control method and system, electronic equipment and a storage medium, and the method comprises the steps: obtaining a physiological state vector based on physiological state data; performing empirical mode decomposition on the physiological state vector, and constructing a biological rhythm digital twinborn model; based on the physiological state data, the environmental parameters and the biological rhythm digital twinborn model, using a PC algorithm of a Pearl causal inference framework to discover a causal structure between the environmental parameters and the physiological state data, and constructing an environmental physiological causal relationship model; constructing a multi-objective optimization problem according to the biological rhythm digital twinborn model and the environmental physiological causal relationship model; the multi-objective optimization problem is converted into a quadratic unconstrained binary optimization function, and a Pareto optimal solution set is obtained for the quadratic unconstrained binary optimization function according to a quantum optimization algorithm; and regulating and controlling the environment and equipment cluster based on the Pareto optimal solution set and the biological rhythm digital twinborn model.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of environment and device cluster regulation, and particularly relates to an environment and device cluster regulation method and system, an electronic device and a storage medium. BACKGROUND

[0002] With the continuous progress of Internet of Things technology and the rapid development of artificial intelligence technology, smart home systems have gradually become an indispensable means in people's daily life. It greatly improves people's quality of life and makes the home environment more intelligent, convenient and comfortable. Through the interconnection of Internet of Things and intelligent analysis of artificial intelligence, the smart home system can realize remote control, automatic management and personalized service of home devices, greatly improving people's life experience and happiness.

[0003] The existing smart home system mainly adopts a control mode based on rules or simple machine learning, and only controls according to environmental parameters or user operation behavior, ignoring the actual needs of human physiological health, and having poor self-adaptation ability. SUMMARY

[0004] The embodiments of the present application provide an environment and device cluster regulation method and system, an electronic device and a storage medium to solve the problems of related technologies. The technical solutions are as follows:

[0005] In a first aspect, the embodiments of the present application provide an environment and device cluster regulation method, comprising:

[0006] Obtaining physiological state data and environmental parameters, and obtaining a physiological state vector based on the physiological state data;

[0007] Performing empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model;

[0008] Discovering a causal structure between the environmental parameters and the physiological state data based on the physiological state data, the environmental parameters and the biological rhythm digital twin model using a PC algorithm of a Pearl causal inference framework to construct an environmental physiological causal relationship model;

[0009] Constructing a multi-objective optimization problem according to the biological rhythm digital twin model and the environmental physiological causal relationship model;

[0010] Converting the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and obtaining a Pareto optimal solution set by optimizing the quadratic unconstrained binary optimization function according to a quantum optimization algorithm;

[0011] Regulating the environment and device cluster based on the Pareto optimal solution set and the biological rhythm digital twin model.

[0012] In a second aspect, the embodiments of the present application provide a regulation system of an environment and a device cluster, comprising:

[0013] A first obtaining module is configured to obtain physiological state data and environment parameters, and obtain a physiological state vector based on the physiological state data;

[0014] A first constructing module is configured to perform empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model;

[0015] A second constructing module is configured to discover a causal structure between the environment parameters and the physiological state data based on the physiological state data, the environment parameters and the biological rhythm digital twin model by using a PC algorithm of a Pearl causal inference framework, and construct an environment-physiology causal relationship model;

[0016] A third constructing module is configured to construct a multi-objective optimization problem according to the biological rhythm digital twin model and the environment-physiology causal relationship model;

[0017] A first obtaining module is configured to convert the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and obtain a Pareto optimal solution set by performing quantum optimization algorithm on the quadratic unconstrained binary optimization function;

[0018] A first regulation module is configured to regulate the environment and the device cluster based on the Pareto optimal solution set and the biological rhythm digital twin model.

[0019] In a third aspect, the embodiments of the present application provide an electronic device, comprising:

[0020] at least one processor; and a memory connected with the at least one processor in communication; wherein

[0021] The memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the method in any one of the embodiments of the above aspects.

[0022] In a fourth aspect, the embodiments of the present application provide a computer-readable storage medium, which stores computer instructions, and when the computer instructions are run on a computer, the method in any one of the embodiments of the above aspects is performed.

[0023] The advantages or beneficial effects of the above technical solutions at least include:

[0024] In the embodiment, the method for regulating the environment and device cluster achieves accurate prediction of the user's future physiological state by constructing a personalized biological rhythm digital twin model, enabling the environment and device cluster to adjust the environment in advance, from passive response to active prediction; the PC algorithm of the Pearl causal inference framework is used to discover the causal structure between the environmental parameters and the physiological state data, and an environmental physiological causal relationship model is constructed to provide a scientific basis for environmental regulation, avoiding blind regulation based on experience or correlation; and quantum optimization technology is used to achieve global optimal balance of multiple objectives, achieving optimization in multiple dimensions such as health, comfort, and energy saving, breaking through the limitations of traditional optimization algorithms; realizing the technical leap from passive response to active prediction, from general control to personalized adjustment, and from single optimization to multi-objective balance. It can automatically adjust the environment and device cluster according to the user's environmental parameters and physiological state data in different scenarios, to adapt to the user's physiological health needs, follow the adaptive biological rhythm, predict in advance, adjust in real time, and have strong self-adaptation and self-regulation capabilities, thereby solving the technical problems of existing smart home systems that mainly use rule-based or simple machine learning control methods, and only control based on environmental parameters or user operation behavior, ignoring the actual needs of human physiological health, and having poor self-adaptation capabilities.

[0025] The above summary is intended to illustrate only and is not intended to be limiting in any way. Further aspects, implementations, and features of the present application will be apparent to those skilled in the art from the description of the illustrative aspects, implementations, and features described above in connection with the accompanying drawings and the following detailed description. BRIEF DESCRIPTION OF DRAWINGS

[0026] In the drawings, like reference numerals refer to like elements throughout the various drawings. The drawings are not necessarily to scale, the emphasis instead being placed upon illustrating certain principles of the application. It should be understood that the drawings are merely depictions of some embodiments of the application and should not be construed as limiting the scope of the application.

[0027] Figure 1 A flowchart of a method for regulating an environment and device cluster according to an embodiment of the present application.

[0028] Figure 2 A block diagram of an electronic device according to an embodiment of the present application. DETAILED DESCRIPTION

[0029] In the following, only some exemplary embodiments are described in brief. As will be appreciated by those skilled in the art, the described embodiments can be modified in various different ways without departing from the spirit or scope of the present application. Therefore, the drawings and the description are to be considered as merely illustrative and not restrictive.

[0030] Figure 1 A flow chart of an environment and device cluster regulation method according to an embodiment of the present application is shown. As shown in the figure, an environment and device cluster regulation method can include: Figure 1

[0031] S110: Obtain physiological state data and environment parameters, and obtain a physiological state vector based on the physiological state data;

[0032] S120: Perform empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model;

[0033] S130: Discover the causal structure between the environment parameters and the physiological state data based on the physiological state data, the environment parameters and the biological rhythm digital twin model using the PC algorithm of the Pearl causal inference framework, and construct an environment-physiology causal relationship model;

[0034] S140: Construct a multi-objective optimization problem according to the biological rhythm digital twin model and the environment-physiology causal relationship model;

[0035] S150: Convert the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and obtain a Pareto optimal solution set by optimizing the quadratic unconstrained binary optimization function according to a quantum optimization algorithm;

[0036] S160: Regulate the environment and device cluster based on the Pareto optimal solution set and the biological rhythm digital twin model.

[0037] ​In the embodiment, the method for regulating the environment and device cluster achieves accurate prediction of the user's future physiological state by constructing a personalized biological rhythm digital twin model, enabling the environment and device cluster to adjust the environment in advance and changing from passive response to active prediction; the PC algorithm of the Pearl causal inference framework is used to discover the causal structure between the environmental parameters and the physiological state data, and an environmental physiological causal relationship model is constructed to provide a scientific basis for environmental regulation, avoiding blind regulation based on experience or correlation; and quantum optimization technology is used to achieve global optimal balance of multiple objectives, achieving optimization in multiple dimensions such as health, comfort, and energy saving, breaking through the limitations of traditional optimization algorithms; realizing the technical leap from passive response to active prediction, from general control to personalized regulation, and from single optimization to multi-objective balance. The method can automatically adjust the environment and device cluster according to the user's environmental parameters and physiological state data in different scenarios, to adapt to the user's physiological health needs, and has strong capabilities of early prediction, real-time adjustment, self-adaptation, and self-regulation, thereby solving the technical problem that existing smart home systems mainly use rule-based or simple machine learning control methods, and only control according to environmental parameters or user operation behavior, ignoring the actual needs of human physiological health, and having poor self-adaptation capability.

[0038] The method for regulating the environment and device cluster of the embodiment can be applied in the scenario of a home environment, and can also be used in office places such as companies, conference rooms, libraries, or study rooms, and can also be adapted for use in hotels and restaurants, and has a wide range of applications.

[0039] In step S110, physiological state data and environmental parameters are obtained, and a physiological state vector is obtained based on the physiological state data.

[0040] In the embodiment, the physiological state data includes heart rate HR(t), heart rate variability HRV(t), body temperature BT(t), systolic blood pressure SBP(t), and diastolic blood pressure DBP(t), and can be characterized by a comprehensive evaluation vector of autonomic nervous state, a body temperature rhythm parameter vector, and a sleep state vector. The physiological state vector can be obtained by integrating the comprehensive evaluation vector of autonomic nervous state, the body temperature rhythm parameter vector, and the sleep state vector.

[0041] Based on the piezoelectric sensor array, the ballistocardiogram signal is collected and the heart rate variability is analyzed. When the heart beats, the mechanical vibration generated by the blood pump will be transmitted to the body surface through the body tissue. When a person lies on a bed, this vibration will be transmitted to the mattress, causing a small pressure change. This phenomenon is called ballistocardiogram (BCG).

[0042] The present application captures this weak signal through a piezoelectric thin film sensor array. Piezoelectric materials have the property of converting mechanical energy to electrical energy and vice versa - when subjected to pressure, they generate an electric charge, the amount of which is proportional to the pressure.

[0043] Specifically, the piezoelectric thin film sensor array arrangement employs a matrix structure of M x N piezoelectric thin film sensors, for example: an 8 x 6 array. This design is based on the following considerations:

[0044] When a person lies down, the main pressure is distributed in the head, shoulder, hip and leg areas. For example: an 8 x 6 array arranged with a 15 cm spacing can cover an area of 120 cm x 90 cm, which is sufficient to cover the lying range of most adults. The 15 cm spacing is the result of optimization. Too small a spacing will increase the cost and computational complexity, and too large a spacing may miss critical signals. This spacing ensures that at least 3-4 sensors can capture the signals in the heart area. The 48 sensors provide sufficient redundancy. Even if some sensors fail or have poor signal quality, the system can still obtain sufficient information through other sensors.

[0045] M x N piezoelectric thin film sensors are evenly arranged under the mattress to form a sensor array P:

[0046] P = {p ij |i∈[1,M],j∈[1,N]}

[0047] Where: M = 8 (number of longitudinal sensors), N = 6 (number of transverse sensors), sensor spacing: d = 15 cm, sampling frequency: fs = 1000 Hz;

[0048] The original pressure signal p ij collected by the (i, j) position sensor at time t is:

[0049] position(p ij ) = (d x (i - 1), d x (j - 1)) cm

[0050] When the heart contracts, blood is pumped into the aorta, generating a counterforce that causes a slight displacement of the body. This displacement, although invisible to the naked eye (usually less than 1 mm), is sufficient to be detected by sensitive piezoelectric sensors. The composition of the original signal can be represented as:

[0051] V ij (t) = V BCG (t) + V resp (t) + V motion (t) + V noise (t)

[0052] Where: V BCG(t): Ballistocardiogram component, which is the target signal we want to extract, V resp (t): Respiratory motion induced signal, lower frequency (0.2-0.5 Hz), V motion (t): Body motion induced signal, usually irregular large amplitude variation, V noise (t): Environmental and electronic noise.

[0053] The ballistocardiogram component can be further modeled as:

[0054]

[0055] where: N beats : Number of heartbeats within the time window, A k : Amplitude of the kth heartbeat, t k : Time of the kth heartbeat, h(t): Waveform template function of a single heartbeat.

[0056] Butterworth filter is chosen because it has the flattest frequency response within the passband and does not introduce signal distortion. A fourth order design provides a good balance of roll-off characteristics and computational efficiency. The transfer function of the filter is derived in detail: For a fourth order Butterworth bandpass filter, it can be decomposed into a cascade of a second order lowpass and a second order highpass:

[0057] Lowpass part (cutoff frequency fH= 40 Hz):

[0058]

[0059] Highpass part (cutoff frequency f L = 0.5 Hz):

[0060]

[0061] Total transfer function:

[0062]

[0063] where:

[0064] ω L = 2πf L = 3.14 rad / s, ω H = 2πf H = 251.3 rad / s

[0065] The filtering process is represented in time domain as a convolution operation:

[0066]

[0067] where h(t) is the filter impulse response, which is obtained by inverse Laplace transform of H(s).

[0068] R-wave is the most prominent peak in ECG, corresponding to ventricular depolarization. In BCG signal, R-wave corresponds to the maximum pressure change. The core idea of adaptive threshold algorithm is: the threshold is not fixed, but dynamically adjusted according to the local statistical properties of the signal. This solves the problem of signal amplitude changing over time (e.g. user changing posture). The calculation process of dynamic threshold includes:

[0069] Sliding window statistics calculation: local mean μ w (t) :

[0070]

[0071] Local standard deviation σ w (t) :

[0072]

[0073] Where W = 2000 sample points (corresponding to 2 seconds), the choice of this window length is based on: need to contain at least 2-3 heartbeat periods to obtain stable statistics, but cannot be too long, otherwise it cannot adapt to the rapid changes of the signal.

[0074] Threshold setting:

[0075] Th(t) = μ w (t) + k · σ w (t)

[0076] The coefficient k = 0.6 is the result of experimental optimization, balancing the detection sensitivity and false positive rate.

[0077] Peak detection and verification, candidate R-wave detection:

[0078] R candidate = {t | V filtered (t) > Th(t)}

[0079] Local maximum verification:

[0080]

[0081] Where Δ = 50 milliseconds, to ensure that the detected is a real peak and not a noise spike.

[0082] For heart rate variability (HRV), it is not only the statistical quantity of heartbeat interval, it reflects the balance state of the autonomic nervous system:

[0083] The meaning of SDNN (standard deviation):

[0084]

[0085] SDNN reflects the overall activity of the autonomic nervous system. The higher the value, the stronger the autonomic regulation, the better the body's adaptability. The SDNN of normal adults is usually in the range of 30-60 milliseconds.

[0086] Significance of RMSSD (Root Mean Square of Successive Differences):

[0087]

[0088] RMSSD mainly reflects the activity of the parasympathetic nervous system (vagus nerve). The parasympathetic nervous system is responsible for "rest and digestion" functions, and an increase in RMSSD indicates that the body is in a relaxed state.

[0089] Significance of pNN50 (Percent of Normal-to-Normal intervals greater than 50ms):

[0090]

[0091] pNN50 also reflects parasympathetic nervous activity, but is more sensitive to rapid changes.

[0092] Physiological basis of frequency domain indicators: convert RR interval series to frequency domain by Fast Fourier Transform (FFT):

[0093]

[0094] Power spectral density calculation uses the Welch method to reduce noise:

[0095]

[0096] where K is the number of segments, each with 50% overlap. Low frequency power (LF, 0.04-0.15 Hz):

[0097]

[0098] LF power reflects the combined action of sympathetic and parasympathetic nerves, mainly related to blood pressure regulation. High frequency power (HF, 0.15-0.4 Hz):

[0099]

[0100] HF power purely reflects parasympathetic nervous activity, synchronized with respiratory rhythm, and is an indicator of autonomic nervous balance:

[0101]

[0102] This ratio reflects the balance between sympathetic and parasympathetic nerves. An increase in the ratio indicates that the sympathetic nervous system is dominant (stress state), and a decrease in the ratio indicates that the parasympathetic nervous system is dominant (relaxed state).

[0103] The final output is a comprehensive assessment vector of the user's autonomic nervous state HRV(t):

[0104] HRV(t) = [SDNN(t), RMSSD(t), pNN50(t), LF / HF(t)] T

[0105] Through the comprehensive assessment vector of the user's autonomic nervous state HRV(t), the physiological indicators of the user's autonomic sympathetic nervous state can be understood, such as the understanding of the heart rate situation, and can be applied to the analysis and processing of subsequent steps.

[0106] The human body is a constant temperature system, and the core body temperature is maintained at about 37℃, but there is a circadian rhythm fluctuation. This fluctuation is controlled by the suprachiasmatic nucleus of the hypothalamus (the biological clock center), and the amplitude is about 1-1.5℃. Body temperature rhythm is one of the most reliable circadian rhythm markers. The present application adopts infrared thermal imaging technology, based on Planck's blackbody radiation law, and calculates the body temperature by measuring the infrared radiation emitted by the human body.

[0107] According to the Stefan-Boltzmann law, the total radiation power P of a black body is:

[0108] P = σ·A·T 4

[0109] Where: σ = 5.67 × 10 -8 W / (m 2 ·K 4 ), Stefan-Boltzmann constant, A: radiating surface area, T: absolute temperature.

[0110] The human skin is not an ideal black body, so the emissivity correction needs to be introduced, and the total radiation power P skin of the corrected black body is:

[0111]

[0112] Where, the emissivity of human skin ε ≈ 0.98, close to black body.

[0113] In addition, the radiation power received by the infrared camera also needs to consider several factors:

[0114] 1. Distance attenuation: according to the inverse square law of radiation intensity:

[0115]

[0116] Where d is the distance from the target to the camera.

[0117] 2. Atmospheric absorption:

[0118] Infrared radiation is absorbed by water vapor and carbon dioxide when propagating in the atmosphere:

[0119] τ atm = e -α·d

[0120] In the 8-14 pm long-wave infrared window, the absorption coefficient a ~ 0.006 m -1 .

[0121] 3. Lens collection efficiency:

[0122]

[0123] Where A lens is the effective area of the lens.

[0124] Taking these factors into account, the power received by the camera is:

[0125]

[0126] In the actual implementation of this embodiment, the human body temperature needs to be obtained, and the temperature is specifically inverted from the received power as follows:

[0127] Ambient radiation compensation, the total radiation received by the camera includes target radiation and environmental reflection:

[0128] P total = ε·P target + (1-ε)·P environment

[0129] Measure the ambient temperature T env and compensate:

[0130]

[0131] Multi-point measurement and weighted average, each pixel in the infrared image corresponds to a temperature value. For the face area, a weighted average is used:

[0132]

[0133] The weight w i is set according to the pixel position, and the weight of the forehead area is the highest (the blood vessels are rich, and the correlation with the core body temperature is the strongest).

[0134] The relationship between skin temperature and core body temperature is affected by multiple factors, including:

[0135] Linear regression model:

[0136] T core (t) = a0 + a1·T forehead (t) + a2·T ambient (t) + a3·RH(t) + a4·v air (t)

[0137] Wherein: T forehead (t): Forehead temperature, obtained through infrared measurement, T ambient (t): Ambient temperature, obtained through an environmental sensor; RH(t): Relative humidity, affecting skin evaporative heat dissipation; v air (t): Air velocity, which affects convective heat dissipation.

[0138] The regression coefficients were obtained through training with large-scale clinical data:

[0139] a h =25.43, intercept term; a1=0.384, forehead temperature coefficient; a2=-0.042, ambient temperature correction coefficient; a3=-0.018, humidity correction coefficient; a6=-0.025, wind speed correction coefficient. The R-value of this model is... 2 The value reached 0.92, and the average error was less than 0.2℃.

[0140] The body temperature rhythm exhibits an approximate cosine waveform. Using the least squares method for fitting, the objective function is:

[0141]

[0142] This is a nonlinear least squares problem, solved using the Levenberg-Marquardt algorithm:

[0143] Initial value estimation:

[0144]

[0145] After iterative update:

[0146] θ (k+1) =θ (k) -(J T J+λI) -1 J T r

[0147] in:

[0148] θ=[T mesor A temp ,φ temp ] T : Parameter vector; J: Jacobian matrix; rr: Residual vector; λ: Damping parameter.

[0149] Convergence criteria:

[0150] ||θ (k+1) -θ (k) ||<10 -6

[0151] The extracted rhythm parameters have important physiological significance:

[0152] T mesor Mean body temperature reflects basal metabolic rate; A temp Amplitude, reflecting rhythm intensity (approximately 0.5℃ in healthy individuals); φ temp Phase: Reflects the phase of the biological clock (normal human body temperature is lowest between 4-6 am).

[0153] The final output is the body temperature rhythm parameter vector TEMP(t):

[0154] TEMP(t)=[T mesor (t), A temp (t), φ temp (t)] T

[0155] The body temperature rhythm parameter vector TEMP(t) can characterize the user's body temperature rhythm.

[0156] Sleep is not a uniform state, but a periodic process consisting of multiple stages:

[0157] Wake period: EEG is dominated by beta waves (13-30Hz);

[0158] Light sleep (N1-N2): Theta waves (4-8Hz) increase, and sleep spindle waves appear;

[0159] Deep sleep (N3): Delta waves (0.5-4Hz) dominate, growth hormone secretion;

[0160] Rapid eye movement (REM) sleep: brain waves similar to those during wakefulness, the dreaming period.

[0161] The breathing pattern is different in each sleep stage:

[0162] Awake: Irregular, easily controlled by consciousness;

[0163] Light sleep: regular, with a slightly lower frequency;

[0164] Deep sleep: deep and slow, very regular;

[0165] REM sleep: Irregular, with occasional apnea.

[0166] This application infers sleep stages by analyzing the characteristics of breathing sounds, avoiding the complexity of traditional polysomnography. Specifically, it infers sleep stages by collecting breathing sounds through a four-element microphone array.

[0167] Configure a four-element microphone array, determine spatial localization capabilities, and locate the sound source using the time difference of arrival (TDOA):

[0168]

[0169] in:

[0170] x s : Location of the sound source, x i ,x j : Microphone position, c = 343 m / s: Speed ​​of sound.

[0171] Noise is suppressed, and beamforming technology is used to enhance the target direction signal.

[0172]

[0173] Where the weight w i and delay τ i Optimize based on the direction of the sound source.

[0174] By configuring four microphones to form a quad microphone array, it can provide multi-angle sound acquisition and redundancy, so that even if one fails, it can still maintain basic functions.

[0175] Breath sound feature extraction can be performed using Short Time Fourier Transform (STFT) for time-frequency analysis:

[0176]

[0177] Where: N = 256: window length (5.3ms@48kHz); H = 128: step length (50% overlap); w[n]: Hamming window function.

[0178] The Hamming window was chosen because it provides a good balance between main lobe width and side lobe suppression.

[0179]

[0180] The calculation of respiratory band energy includes:

[0181]

[0182] Where k1 and k2 correspond to frequency indices of 0.2 to 0.5 Hz.

[0183] Calculation of respiratory energy variance:

[0184]

[0185] Respiratory Energy Variance (VAR) breath (t) reflects the regularity of breathing, with the smallest variance during deep sleep.

[0186] Zero-crossing rate ZCR(t) calculation:

[0187]

[0188] The zero-crossing rate ZCR(t) reflects the high-frequency components of the signal; ZCR decreases during snoring.

[0189] Calculation of spectral centroid SC(t):

[0190]

[0191] The centroid SC(t) reflects the "center of gravity" of the spectrum, and the centroid decreases as breathing depth increases.

[0192] Mel frequency cepstral coefficients (MFCC) are used to capture vocal tract features and distinguish between normal breathing and abnormal events (such as sleep apnea).

[0193] Eigenvectors of Mel-frequency cepstral coefficients (MFCC):

[0194] F sleep (t)=[E breath VAR breath ,ZCR,SC,MFCC1,...,MFC 13 ] T

[0195] Hidden Markov Model (HMM) Sleep Stages

[0196] Hidden Markov Models (HMMs) are particularly well-suited for sleep stages because sleep stages have Markovian properties (the next stage largely depends on the current stage). We observe acoustic features rather than the sleep stages themselves (latent states).

[0197] HMM parameters include:

[0198] State transition probability matrix:

[0199] A = [a ij ] 4×4 ,a ij =P(S) t+1 =j|S t =i)

[0200] Based on the initial understanding of sleep physiology: sleep stages usually progress sequentially, with deep sleep mainly occurring in the first half of the night and REM sleep increasing in the second half of the night.

[0201] Observing the probability distribution, we assume a Gaussian mixture model (GMM):

[0202]

[0203] M = 3 Gaussian components is usually sufficient.

[0204] Initial state distribution:

[0205]

[0206] It is usually assumed that the initial state is conscious: π wake =1.

[0207] The Viterbi algorithm is used to infer the most likely sleep sequence, specifically given an observation sequence O = {o1,...,o...}. T Find the most likely sequence of states:

[0208]

[0209] Viterbi algorithm dynamic programming recursion:

[0210] initialization:

[0211] δ1(i)=π i ·b i (o1)

[0212] ψ1(i)=h

[0213] Further deduction:

[0214]

[0215] termination

[0216]

[0217] Backtracking:

[0218]

[0219] Calculate sleep quality indicators from sleep stage results:

[0220] Sleep efficiency:

[0221]

[0222] TST represents total sleep time, and TIB represents time spent in bed.

[0223] Deep sleep percentage:

[0224]

[0225] For healthy adults, the percentage of deep sleep (N3%) should reach 15-20%.

[0226] REM sleep percentage:

[0227]

[0228] The normal range for REM sleep percentage is 20-25%.

[0229] Sleep fragmentation index:

[0230]

[0231] The sleep fragmentation index reflects sleep continuity; the lower the value, the better.

[0232] The final output is a comprehensive assessment of the sleep state, specifically the sleep state vector SLEEP(t):

[0233] SLEEP(t)=[S(t),SE(t),N3%(t),REM%(t),FI(t)] T

[0234] The sleep state vector SLEEP(t) can characterize the user's sleep state.

[0235] Millimeter-wave radar can simultaneously detect heartbeats and peripheral pulse waves. Blood pressure can be estimated by measuring the time difference between the pulse wave's propagation from the heart to the periphery (such as the carotid artery), i.e., the pulse wave travel time (PTT). Since blood pressure and PTT are inversely proportional:

[0236] Blood pressure estimation formula: SBP = α / PTT + β

[0237] α and β are individualized parameters that need to be obtained through initial calibration.

[0238] In this embodiment, environmental parameters include temperature (TEMP(t), humidity (HUM(t), illuminance (LUX(t), color temperature (CCT(t)), carbon dioxide concentration (CO2(t)), and noise level (NOISE(t)). Specifically, each environmental parameter affects health through a specific physiological pathway.

[0239] Temperature (TEMP(t)): Affects the thermoregulatory center, alters peripheral vasomotor activity, and influences metabolic rate.

[0240] Humidity (HUM(t)): affects sweat evaporation, alters the state of the respiratory mucosa, and affects the survival of pathogens.

[0241] Illuminance LUX(t): inhibits / promotes melatonin secretion, regulates circadian rhythms, and affects mood (seasonal affective disorder).

[0242] Color temperature (CCT(t)) affects melatonin secretion.

[0243] Carbon dioxide concentration (CO2(t)): affects blood oxygen saturation, causes inflammatory responses, and affects cognitive function.

[0244] Noise level NOISE(t): activates the sympathetic nervous system, disrupts sleep, and induces stress response.

[0245] Environment state vector:

[0246] ENV(t) = [T a (t), RH(t), E v (t), CCT(t), CO2(t), PM 2.5 (t),L p (t)] T

[0247] in:

[0248] T a (t): Air temperature (°C); RH(t): Relative humidity (%); E v (t): Illuminance (lux); CCT(t): Color temperature (K); CO2(t): Carbon dioxide concentration (ppm); PM2.5(t): PM2.5 concentration (μg / m³) 3 Lp(t): sound pressure level (dB).

[0249] In step S120, empirical mode decomposition is performed on the physiological state vector to construct a digital twin model of biological rhythms.

[0250] In the embodiments of this application, Empirical Mode Decomposition (EMD) is an adaptive signal decomposition method, particularly suitable for nonlinear and non-stationary physiological signals. Unlike Fourier Transform, EMD does not require preset basis functions, but decomposes the signal based on its own characteristics.

[0251] Physiological state vectors can contain rhythms across multiple time scales, for example:

[0252] Ultra-short rhythms (5–20 minutes): reflect immediate physiological regulation, such as stress response;

[0253] Super-solar rhythms (90–120 minutes): such as sleep cycles and hormonal pulses;

[0254] Circadian rhythm (24 hours): The basic rhythm controlled by the biological clock;

[0255] Long-term trends (>24 hours): Seasons, menstrual cycles, etc.;

[0256] These rhythms overlap, and direct analysis would cause them to interfere with each other.

[0257] In this embodiment, empirical mode decomposition is performed on the physiological state vector, and multi-scale biorhythm decomposition results are obtained by identifying the user's rhythmic characteristics. Then, a coupling network between physiological rhythms is constructed based on the multi-scale biorhythm decomposition results to discover the mutual influence relationships between rhythms. Finally, a physiological state prediction model is established by integrating the multi-scale biorhythm decomposition results and the coupling network between physiological rhythms. There is a clear data dependency between these three sub-steps: the multi-scale biorhythm decomposition results are the foundation, the coupling network between physiological rhythms is the bridge, and the physiological state prediction model is the final goal, ultimately resulting in a biorhythm digital twin model.

[0258] By constructing personalized biorhythm digital twin models based on users' physiological state data, we can accurately predict users' physiological state in the future. This enables the environment and equipment clusters to adjust in advance, shifting from passive response to proactive prediction. This allows for better adaptation to users' physiological state, resulting in a more comfortable environment and a significantly improved user experience.

[0259] In this embodiment, empirical mode decomposition (EMD) is performed on the physiological state vector to obtain multi-scale circadian rhythm decomposition results. By performing EMD on the physiological state vector, the original physiological state data is decomposed into rhythmic components at different time scales. For each identified rhythmic component, it can include the period length of the substructure (e.g., 15 minutes, 24 hours, 7 days), amplitude (reflecting the strength of the rhythm), phase information (reflecting the time shift of the rhythm), and the proportion of that rhythm in the overall physiological variability. For example, for a user with sleep phase delay, their circadian rhythm phase delay might be identified as 3.2 hours, meaning their physiological "biological clock midnight" actually occurs at 3:12 AM. In this case, this phase delay information is crucial for subsequent environmental and device cluster control decisions.

[0260] A coupling network among physiological rhythms was constructed. Granger causality analysis was used to analyze the decomposition results of multi-scale biorhythms. Transfer entropy was used to calculate the decomposition results, and the influence types of the decomposition results were analyzed to construct a coupling network between different physiological indicators. The physiological rhythm coupling network contains multiple nodes (representing different physiological indicators) and directed edges (representing causal relationships). Each edge is labeled with three key attributes: causal strength (Granger causality index, ranging from 0 to 1), time delay (the time required for influence transmission, in minutes), and information transmission amount (transfer entropy value, reflecting the strength of information flow). For example, the physiological rhythm coupling network might find that the causal influence strength of body temperature rhythm on heart rate variability is 0.42, with a delay time of 10 minutes. This indicates that changes in body temperature will affect heart rate variability after 10 minutes, with a moderate influence strength.

[0261] Based on the results of multi-scale rhythm decomposition, the physiological rhythm coupling network, and the physiological state vector, a physiological state prediction model is obtained. This model can include complete physiological state predictions for up to 32 time points (one prediction point every 15 minutes). The prediction for each time point can include: heart rate value, heart rate variability indices (SDNN and RMSSD), core body temperature, sleep stage (wakefulness, N1 light sleep, N2 moderate sleep, N3 deep sleep, REM rapid eye movement sleep), alertness level (a continuous value from 0 to 1), and prediction confidence, among other future data. In addition to conventional time-series predictions, the physiological state prediction model can also identify key physiological events, such as the expected sleep onset time, the peak of deep sleep, and the REM sleep window, and provide probability estimates for these events. This allows for predictions of the user's future physiological state, enabling adaptive adjustments by the environment and equipment clusters to better meet user needs and improve user satisfaction.

[0262] In one example, Empirical Mode Decomposition (EMD) was performed on a user's physiological state vector over seven consecutive days, revealing rhythmic patterns across multiple time scales. Taking the SDNN metric of heart rate variability as an example, EMD decomposes it into 4–8 intrinsic modes (IMFs). In the example below, the physiological state vector is decomposed into 6 intrinsic modes. In other instances, it can also be decomposed into 8 intrinsic modes, depending on the specific circumstances.

[0263] The first intrinsic mode, IMF1 – ultrashort rhythm (5–15 minutes): This component reflects the user's immediate stress response. Analysis revealed two peaks daily, between 2:00–3:00 PM and 7:00–8:00 PM, corresponding to afternoon drowsiness and the post-dinner digestive load period, respectively. The amplitude is approximately 8 ms, accounting for 15% of the total variance. This rapid fluctuation is highly correlated with work-related stress events: when receiving an urgent task, IMF1 rises by 10–15 ms within 5 minutes.

[0264] The second intrinsic modality, IMF2-Supersolar Rhythm (90 minutes), reveals the basal rest-activity cycle (BRAC). Users' average BRAC cycle is 95 minutes, slightly longer than the standard 90 minutes. This rhythm manifests as periodic fluctuations in attention during work: a significant decrease in focus after every 1.5 hours of work. Interestingly, this rhythm translates into sleep cycles during sleep, with the complete N1-N2-N3-REM cycle also approximately 95 minutes.

[0265] The third intrinsic mode, IMF3-semi-diurnal rhythm (12 hours), exhibited an unexpected 12-hour rhythm, peaking at 6 AM and 6 PM. This may be related to the user's bipolar lifestyle: daytime work and late-night leisure activities create two activity peaks. This abnormal bimodal pattern exacerbates circadian rhythm disruption.

[0266] The fourth intrinsic mode, IMF4-circadian rhythm (24 hours): This is the most dominant rhythmic component, accounting for 45% of the total variance. However, phase analysis shows that the user's circadian rhythm phase is delayed by 3 hours. In normal individuals, the SDNN is lowest between 3-4 AM and highest between 3-4 PM; while for this user, the lowest point is delayed until 6-7 AM, and the highest point is between 6-7 PM. This phase delay is the root cause of late nights.

[0267] The fifth intrinsic mode, IMF5-weekly rhythm (2-3 days), reflects the difference between weekdays and weekends. Physiological stress gradually accumulates from Monday to Wednesday; it peaks on Thursday; recovery begins on Friday; and overcompensation (revenge sleep deprivation) occurs on the weekend. This pattern creates a vicious cycle.

[0268] The sixth intrinsic modality, IMF6, shows a long-term trend: In the first month of monitoring, the overall SDNN showed a downward trend (-0.5ms / day), indicating a deterioration in autonomic nervous function. This is an early warning signal that sub-health is progressing towards disease.

[0269] Granger causality analysis was used to determine the results of Granger causality analysis between users' physiological systems.

[0270] Causal relationship between body temperature and heart rate: Granger causality index (GC) = 0.42 (p < 0.001), with a time lag of approximately 10 minutes. Specifically, for every 0.1°C increase in body temperature, heart rate increases by 2-3 bpm. This relationship is weakest at night (GC = 0.25) and strongest during the day (GC = 0.58), reflecting the influence of diurnal metabolic activity.

[0271] The causal relationship between respiration and heart rate variability (HRV) was GC = 0.68 (p < 0.001), indicating near real-time correlation. Deep breathing increased RMSSD, a manifestation of respiratory sinus arrhythmia. Users experienced shallow breathing under stress, leading to a decrease in HRV, creating a vicious cycle. Once the system identified this critical point, environmental interventions could break this cycle.

[0272] Causal relationship between sleep stage and body temperature: GC = 0.35 (p < 0.01), with a time lag of approximately 30 minutes. Body temperature drops by 0.2-0.3℃ after entering deep sleep, which is a normal physiological phenomenon. However, if the user's deep sleep duration is short, the temperature drop may be insufficient, affecting sleep quality.

[0273] The causal relationship between heart rate variability (HRV) and sleep quality: GC = 0.51 (p < 0.001), affecting the following day. HRV levels 2 hours before bedtime can predict sleep quality for the night. When the LF / HF ratio > 3 (sympathetic overactivation), the proportion of deep sleep decreases by 50%. This finding provides a basis for regulating the pre-sleep environment.

[0274] The transfer entropy results were calculated from the decomposition of multi-scale biological rhythms, revealing nonlinear relationships. For example, there is a nonlinear relationship between the phase of the body temperature rhythm and sleep duration: when the lowest point of body temperature is delayed by more than 2 hours, the time to fall asleep increases exponentially. This explains why users find it harder to fall asleep the more they stay up late.

[0275] A personalized LSTM prediction model was trained based on historical physiological data and historical environmental parameters in the physiological state vector, resulting in a physiological state prediction model. The physiological state prediction model takes into account the user's physiological data and environmental parameters from the past 24 hours (96 15-minute time points) and predicts the physiological state for the next 8 hours (32 time points).

[0276] The LSTM prediction model was trained using historical physiological data and historical environmental parameters from the user's 30-day physiological state vector, totaling 43,200 samples. Key findings during training:

[0277] Attention weighting analysis: The physiological state prediction model has learned to "pay attention" to key time points. For predicting the state at 11 PM, the highest attention weights were: 0.15 for coffee intake at 3 PM (detected by a sudden increase in heart rate), 0.12 for dinner at 6 PM (due to an increase in body temperature), and 0.18 for the previous night's sleep quality. This indicates that the model has captured physiological correlations across time scales.

[0278] Prediction accuracy verification: On the test set, the physiological state prediction model's prediction accuracy for users was:

[0279] Heart rate: MAE = 3.2 bpm, RMSE = 4.1 bpm

[0280] Body temperature: MAE=0.15℃, RMSE=0.19℃

[0281] Sleep stage: Accuracy 82%, Cohen's skappa = 0.76

[0282] HRV index: correlation coefficient r = 0.85

[0283] Anomaly pattern recognition: The physiological state prediction model learns to identify abnormal patterns in users. For example, when it detects caffeine intake after 9 PM (which simultaneously increases heart rate and body temperature), the physiological state prediction model predicts a 2.5-hour delay in sleep onset and a 35% reduction in deep sleep.

[0284] In step S130, based on physiological state data, environmental parameters, and a digital twin model of biological rhythms, the PC algorithm of the Pearl causal inference framework is used to discover the causal structure between environmental parameters and physiological state data, and to construct an environmental-physiological causal relationship model.

[0285] In this embodiment, the PC (Peter-Clark) algorithm is the core algorithm of the Pearl causal inference framework, used to discover causal structures from observational data. The PC algorithm is a causal discovery algorithm based on conditional independence tests. Based on the causal structures learned by the PC algorithm within the Pearl causal inference framework, a quantitative model of environmental-physiological causal relationships is established.

[0286] Specifically, based on the multi-scale biological rhythm decomposition results of the physiological state vector, which may include physiological state data such as heart rate HR(t), heart rate variability HRV(t), core body temperature BT(t), systolic blood pressure SBP(t), and diastolic blood pressure DBP(t), and environmental parameters such as temperature TEMP(t), humidity HUM(t), illuminance LUX(t), color temperature CCT(t), carbon dioxide concentration CO2(t), and noise level NOISE(t), and the four rhythm components IMF1 and IMF2 containing the highest frequency components (5-20 minute cycle), IMF3 and IMF4 containing the intermediate frequency components (1-4 hour cycle), IMF5 and IMF6 containing the low frequency components (12-24 hour cycle), and IMF7 and IMF8 containing the lowest frequency components (3-7 day cycle), a variable matrix is ​​constructed, consisting of 5 physiological state data, 6 environmental parameters, and 4 rhythm components. The initial assumption is that there may be causal relationships between all variables, so there are edges between all variables. Based on the variable matrix, an initial undirected graph is constructed, which is a 26×26 initial undirected graph containing 325 possible edges (26×25 / 2).

[0287] For any two variables X and Y (which can be environmental or physiological variables), and a condition set S (a subset of other variables), a physiological rhythm coupling network is used to test whether X and Y are conditionally independent given S. The results of the conditional independence test are obtained, and edges without causal relationships are removed based on these results, resulting in a skeleton graph. V-structures (also called convergent structures) are triples of the form A→C←B, where A and B are not directly connected. All triples (A,C,B) are found where AC and BC are connected, but A and B are not. For each such triple, it is checked whether A and B become dependent (rather than independent) given C. After identifying V-structure triples, the direction propagation rule is used to determine the direction of edges in the skeleton graph other than those of the triples.

[0288] For the secondary edges in the skeleton graph whose direction cannot be determined by the PC algorithm, temporal information and domain knowledge are used for orientation to obtain the secondary edges whose direction cannot be determined by the PC algorithm and their corresponding directions. Combining the skeleton graph, edge directions, and secondary edge directions obtained above, a causal structure can be derived.

[0289] After determining the causal structure, it is necessary to quantify the strength of each causal edge. Using structural equation modeling, the causal effect of each edge A→B is estimated. The standardized causal strength is obtained by calculating the causal strength at different time scales using multi-scale circadian rhythm decomposition results.

[0290] The final output is a standardized causal strength complete environment-physiological causal model, including: 1. A causal graph structure G = (V, E), where V is 26 nodes (6 environmental + 20 physiological rhythm variables), and E is a defined set of causal edges; 2. A causal parameter matrix P, storing the causal strength S for each edge e = (i, j) ∈ E. ij (between 0 and 1), causal delay D ij (minutes) and confidence level C ij (Based on statistical significance); 3. Dynamic causal equations describing the dynamics of the entire environment and equipment cluster:

[0291] dY1 / dt=f1(X,Y)+h1(t)-α1(Y1-Y 10 )+σ1W1, (Heart Rate Equation)

[0292] dY2 / dt=f2(X,Y)+h2(t)-α2(Y2-Y 20 )+σ2W2, (Equation for heart rate variability)

[0293] dY3 / dt=f3(X,Y)+h3(t)-α3(Y3-Y 30 )+σ3W3, ​​(body temperature equation)

[0294] dY4 / dt=f4(X,Y)+h4(t)-α4(Y4-Y 40 )+σ4W4, (Symptom pressure equation)

[0295] dY s / dt=f s (X, Y) + h s (t)-α s (Y s -Y s0 )+σ s W s (Diastolic pressure equation)

[0296] Wherein, dY1 / dt: instantaneous rate of change of heart rate HR, in (beats / minute) / minute; dY2 / dt: instantaneous rate of change of heart rate variability HRV, in milliseconds / minute; dY3 / dt: instantaneous rate of change of body temperature BT, in °C / minute; dY4 / dt: instantaneous rate of change of systolic blood pressure SBP, in mmHg / minute; dY5 / dt: instantaneous rate of change of diastolic blood pressure DBP, in mmHg / minute.

[0297] Six environmental variables: X1 = TEMP, X2 = HUM, X3 = LUX, X4 = CCT, X5 = CO2, X6 = NOISE; five physiological variables: Y1 = HR, Y2 = HRV, Y3 = BT, Y4 = SBP, Y5 = DBP.

[0298] h1(t) Heart rate rhythm: Diurnal rhythm: A 11 = 5 times / minute, T 11 = 24.1 hours, φ 11 =15 (peak at 3 PM) Supersolar rhythm: A 12 = 2 times / minute, T 12 = 90 minutes (REM sleep cycle), h1(t) = 5 × sin(2π(t-15) / 24.1) + 2 × sin(2π(t-0) / 1.5);

[0299] h2(t) Heart rate variability rhythm: diurnal rhythm: A 21 =10ms, T 21 =24 hours, φ 21 =3 (peak at 3 AM, parasympathetic dominance), h2(t) = 10 × sin(2π(t-3) / 24);

[0300] h3(t) body temperature rhythm: diurnal rhythm: A 31 =0.5℃, T 31 =24.2 hours (Peak at 4 PM), h3(t) = 0.5 × sin(2π(t-16) / 24.2);

[0301] h4(t) Systolic blood pressure rhythm: Diurnal rhythm: A 41 =10mmHg, T 41 =24 hours (Peak at 2 PM), h4(t) = 10 × sin(2π(t-14) / 24);

[0302] h5(t) diastolic blood pressure rhythm: diurnal rhythm: A 51 =5mmHg, T 51 =24 hours (Peak at 2 PM), h5(t) = 5 × sin(2π(t-14) / 24).

[0303] α1 = 0.05 / minute (heart rate), meaning: the deviation of heart rate from baseline by 5% per minute; Y 10 =70 beats / minute (resting heart rate); α2 = 0.03 / minute (heart rate variability), meaning: HRV recovers slowly, recovering 3% per minute; Y 20= 40ms (healthy HRV level); α3 = 0.02 / minute (body temperature), meaning: body temperature regulation is the slowest, recovering 2% per minute; Y 30 =36.8℃ (normal body temperature); α4 = 0.04 / min (systolic blood pressure), meaning: moderate blood pressure regulation speed; Y 40 =120 mmHg (normal systolic blood pressure); α5 = 0.04 / min (diastolic blood pressure), meaning: a rate of adjustment similar to systolic blood pressure; Y 50 =80 mmHg (normal diastolic blood pressure)

[0304] σi: Noise intensity (standard deviation) W of the i-th physiological variable i (t): Standard white noise process (mean 0, variance 1); Noise intensity of each variable: σ1 = 2 beats / minute (random fluctuation of heart rate); σ2 = 3 ms (random fluctuation of HRV); σ3 = 0.05℃ (random fluctuation of body temperature); σ4 = 3 mmHg (random fluctuation of systolic blood pressure); σ5 = 2 mmHg (random fluctuation of diastolic blood pressure).

[0305] In this embodiment, the Pearl causal framework is combined with dynamical systems theory to discover causal structures and determine time delays. The introduction of rhythm-aware causal analysis, considering the modulating effect of biological rhythms on causal relationships, makes it more adaptable to users with different personalities, exhibiting strong adaptability and higher versatility.

[0306] Key physiological events predicted by the physiological state prediction model (such as expected sleep onset time and deep sleep window) can guide the Pearl causal inference framework's PC algorithm to focus on analyzing environmental-physiological causal relationships during these time periods. For example, if the physiological state prediction model predicts that a user will fall asleep at 1:45 AM, the Pearl causal inference framework's PC algorithm will focus on analyzing which environmental factors significantly affect sleep initiation within the time window of 1:15 AM to 2:15 AM. Through this targeted analysis, the Pearl causal inference framework's PC algorithm can discover strong causal relationships within specific time windows, such as "within 30 minutes before sleep onset, for every 1°C decrease in ambient temperature, the sleep latency shortens by 8 minutes."

[0307] Take the causal chain of temperature-sleep quality as an example:

[0308] Causal path analysis revealed two main paths:

[0309] Temperature → Body temperature regulation → Melatonin secretion → Sleep quality (direct pathway);

[0310] Temperature → Comfort level → Stress level → HRV → Sleep quality (indirect path);

[0311] Through path analysis, the total causal effect is decomposed into:

[0312] Direct effect: -0.35 (deep sleep decreases by 3.5% for every 1°C increase in temperature);

[0313] Indirect effect: -0.18 (through stress response);

[0314] Total effect: -0.53;

[0315] This means that lowering the bedroom temperature from 26°C to 22°C can increase deep sleep by 21.2%. Actual verification shows that the prediction error is less than 15%.

[0316] Cross-convergence mapping (CCM) analysis revealed the time lag in the physiological effects of different environmental factors:

[0317] Temperature → Heart Rate: Time delay 8 ± 2 minutes. This is the time it takes for skin temperature receptors to affect the cardiovascular center via nerve conduction.

[0318] Light exposure → melatonin: Time delay 45 ± 10 minutes. It takes time for light signals to inhibit melatonin synthesis through the retina-pineal gland pathway.

[0319] Humidity → Respiration: Time delay 12±3 minutes. The physiological adaptation time of the respiratory mucosa to changes in humidity.

[0320] CO2 → Cognition: Time delay 20 ± 5 minutes. Changes in blood CO2 concentration affect the timing of cerebral blood flow and oxygen supply.

[0321] These time delays are crucial for predictive control. For example, environmental and device clusters can begin reducing light intensity 45 minutes before a user prepares to sleep, ensuring that melatonin levels are optimal at the time of sleep onset, helping users fall asleep comfortably and improving sleep quality.

[0322] In step S140, a multi-objective optimization problem is constructed based on the biological rhythm digital twin model and the environmental physiological causal relationship model.

[0323] In this embodiment, the multi-objective optimization problem includes: a health promotion objective function, a comfort objective function, an energy consumption cost objective function, and an equipment lifespan objective function.

[0324] The health promotion objective function H(X) is used to measure the degree of improvement in user health caused by the environmental control sequence X. Here, X is the environmental control sequence for the next 24 hours, which includes the set values ​​of 6 environmental parameters for each control period.

[0325] First, a physiological state prediction model is used to predict the physiological state trajectory under environmental control X. For the control sequence X = [X(t1), X(t2), ..., X(t...], ... 24 )], where X(t)i The physiological trajectory Y is predicted based on the six environmental parameter settings for the i-th hour. pred =[Y(t1),Y(t2),...,Y(t)] 24 )).

[0326] Then, define the health deviation. For each time-space physiological state Y(t) = [HR(t), HRV(t), BT(t), SBP(t), DBP(t)], calculate the deviation from the target health state:

[0327] Deviation = Σ j w j ×|Y j (t)-Y j-target | / σ j

[0328] Where: Y j_target It is the target value of the j-th physiological indicator (e.g., HR). target =70), σ j w is the standard deviation of the j-th indicator, used for normalization. j These are weights, reflecting the importance of the indicator. The weights are determined based on medical knowledge and the user's health status. For example, for users at risk of hypertension, the weights for blood pressure (w4=0.3, w5=0.3) are relatively high; for users under excessive stress, the weight for HRV (HRV) (w2=0.4) is relatively high.

[0329] Based on the results of multi-scale biological rhythm decomposition, the health function also needs to consider rhythm synchronicity. Ideal environmental control should enhance the regularity of physiological rhythms.

[0330] Circadian rhythm synchronicity calculation: Extract the diurnal rhythm components of the predicted physiological trajectory, calculate rhythm intensity: rhythm component energy / total energy, calculate phase stability: standard deviation of phase changes over several consecutive days, and the complete form of the health promotion objective function H(X):

[0331] H(X) = -[α1×Σ t Health deviation (t) + α2 × (1 - rhythm synchronicity) + α3 × number of abnormal events

[0332] The negative sign is used here because optimization problems are usually defined as minimization, while we aim to maximize health. α1, α2, and α3 are weighting coefficients, typically α1 = 0.5, α2 = 0.3, and α3 = 0.2. The number of abnormal events refers to the number of times physiological indicators exceed the safe range, indicating whether abnormalities will occur under control X using an environmental-physiological causal relationship model.

[0333] The comfort objective function C(X) is used to evaluate the impact of environmental settings on users' subjective comfort. Comfort encompasses multiple dimensions:

[0334] 1. Thermal comfort: Based on the PMV (Predicted Average Voting) model, PMV = f(temperature, humidity, wind speed, clothing, activity level). For indoor environments, this is simplified to:

[0335] Thermal comfort deviation = |TEMP - TEMP preferred |+0.3×|HUM-HUM preferred |

[0336] Among them TEMP preferred and HUM preferred It is the user's preferred temperature and humidity, learned from historical data.

[0337] 2. Visual comfort: Illumination deviation = |LUX - LUX optimal (t)| / LUX optimal (t)

[0338] LUX optimal (t) represents the optimal illuminance at time t, which varies over time: Morning (6-9 am): 300-500 lux, promoting wakefulness; Daytime (9-6 pm): 500-800 lux, maintaining alertness; Evening (6-10 pm): 200-300 lux, preparing for rest; Nighttime (10-6 am): <50 lux, promoting sleep.

[0339] 3. The impact of environmental change rate on comfort

[0340] Using causal delay information from an environmental-physiological causal model, we consider the impact of the rate of environmental change on comfort. Rapid environmental changes can cause discomfort.

[0341] Rate of change penalty = Σ i Σ i β i ×|X i (t+1)-X i (t)|

[0342] Where, β i The sensitivity coefficients for changes in the i-th environmental parameter are: temperature change sensitivity coefficient β1 = 2.0 (rapid temperature changes are the most uncomfortable); illuminance change sensitivity coefficient β3 = 1.5 (sudden changes in light are glaring); humidity change sensitivity coefficient β2 = 0.5 (humidity changes are less perceptible).

[0343] The complete form of the comfort objective function C(X):

[0344] C(X) = γ1 × thermal comfort deviation + γ2 × illumination deviation + γ3 × rate of change penalty

[0345] The weights γ1 = 0.4, γ2 = 0.3, and γ3 = 0.3 can be adjusted based on user feedback.

[0346] The energy consumption objective function E(X) calculates the energy cost required to achieve environmental control X.

[0347] Air conditioning energy consumption model:

[0348] P ac =COP×|TEMP set -TEMP outdoor |×Volume×ρ×c / η

[0349] Where: COP is the coefficient of performance for cooling / heating (typical value 3.5); TEMP outdoor It is the outdoor temperature (obtained from the weather forecast); Volume is the room volume (e.g., 50m²). 3 ); ρ is the air density (1.2 kg / m³). 3 c is the specific heat of air (1000 J / kg·K); η is the equipment efficiency (0.8).

[0350] fficacy

[0351] P light =LUX × Area / luminous efficacy

[0352] Where Area is the room size (e.g., 20m²) 2 ), luminous efficacy It refers to luminous efficacy (100lm / W for LED).

[0353] Humidification / dehumidification energy consumption:

[0354] P hum =|HUM set -HUM current |×Volume×Power perpercent

[0355] Calculate the actual electricity cost by combining time-of-use pricing:

[0356] Electricity cost = Σ t [P total (t)×Price(t)×Δt]

[0357] Where Price(t) is the electricity price at time t:

[0358] Therefore, the energy consumption objective function E(X) is:

[0359] E(X)=Σ t [P ac (t)+Plight (t)+P hum [(t)]×Price(t)

[0360] The equipment life function L(X) is used to evaluate the impact of control strategies on equipment life.

[0361] Start-stop cycle loss: Each start-stop cycle shortens the equipment's lifespan.

[0362] Start-stop loss = Σ i n i ×cost perstarstart

[0363] Where n i It is the number of times device i is started and stopped, cost perstart It is the lifespan cost of each start-up and shutdown.

[0364] Operational stress loss: Equipment aging is accelerated under extreme operating conditions.

[0365] Strength loss = Σ t f(Operating power / Rated power)

[0366] When the operating power approaches the rated power, the f value increases sharply.

[0367] Therefore, the objective function for equipment lifespan is L(X):

[0368] L(X) = δ1 × start-up and shutdown losses + δ2 × intensity losses

[0369] The four objective functions are integrated into a multi-objective optimization problem: minimizeF(X);

[0370] minimize F(X)=[f1(X), f2(X), f3(X), f4(X)]

[0371] Where: f1(X) = -H(X) (health promotion, take the negative sign to minimize); f2(X) = C(X) (comfort deviation); f3(X) = E(X) (energy consumption cost); f4(X) = L(X) (equipment depreciation).

[0372] Constraints:

[0373] Environmental parameter range: 18≤TEMP≤28, 30≤HUM≤70, etc.

[0374] Rate of change limit: |X i (t+1)-X i (t)|≤ΔX imax

[0375] Equipment capacity: P total(t) ≤P max

[0376] The multi-objective optimization problem constructed by the above four objective functions is generally applicable to the control scenarios of environment and equipment clusters. It comprehensively considers the four core dimensions of health, comfort, energy consumption and equipment wear and tear, and can save energy and protect the environment while taking into account the user experience, and ensure the service life of the equipment.

[0377] In step S150, the multi-objective optimization problem is transformed into a quadratic unconstrained binary optimization function, and the Pareto optimal solution set is obtained by applying the quantum optimization algorithm to the quadratic unconstrained binary optimization function.

[0378] In this embodiment, the quantum optimization algorithm processes binary variables, encoding continuous environmental parameters into binary. For a 24-hour control sequence, there are 6 parameters per hour, each encoded with 4 bits, requiring a total of 24 × 6 × 4 = 576 binary variables. This is denoted as a vector q = [q1, q2, ..., q...]. 576 ].

[0379] The multi-objective optimization problem is optimized by performing a second-order unconstrained binary optimization (QUBO) as follows:

[0380] The health promotion objective function, comfort objective function, energy cost objective function, and equipment life objective function in the above multi-objective optimization problem are each transformed into a quadratic unconstrained binary optimization function. For each objective function, its QUBO matrix Q is constructed. i :

[0381] The QUBO (Quadratic Unconstrained Binary Optimization) matrix Q is a symmetric matrix, and its element Qij is defined as follows:

[0382] When i ≠ j:

[0383]

[0384] When i = j:

[0385]

[0386] Among them, f comfort Let f be the comfort objective function. energy Let w1 and w2 be the energy consumption objective function, and w1 and w2 be the weighting coefficients, satisfying w1 + w2 = 1; xi and xj are the decision variables.

[0387] Q1: QUBO matrix of health promotion objective function (576×576); Q2: QUBO matrix of comfort objective function; Q3: QUBO matrix of energy consumption cost objective function; Q4: QUBO matrix of equipment life objective function.

[0388] The computation of matrix elements involves a large number of symbolic operations, but the principle is to expand each term of the original objective function into a product of binary variables, resulting in the QUBO matrix Q in the multi-objective optimization problem. i It is an unconstrained optimization, where the constraints are transformed into penalty terms and added to the QUBO matrix Q. i In this process, we obtain the QUBO matrix q. i It is unconstrained optimization.

[0389] The standard form of the QUBO problem is: minimize f(q) = Σ ij Q ij q i q j

[0390] Where q is a binary vector, each q i It takes the value 0 or 1, and Q is the coefficient matrix.

[0391] Configure QUBO matrix q i The weights of the objective functions for health promotion, comfort, energy cost, and equipment lifespan, for example:

[0392] w1 (health weight): from 0.1 to 0.7, with a step size of 0.1;

[0393] w2 (comfort weight): from 0.1 to 0.7, with a step size of 0.1;

[0394] w3 (energy consumption weight): from 0.1 to 0.7, with a step size of 0.1;

[0395] w4 (device weight): 1-w1-w2-w3 (guaranteed sum to 1);

[0396] For each weight combination (w1, w2, w3, w4):

[0397] Construct a combinatorial QUBO matrix Q combined :

[0398] Q combined = w1×Q1+w2×Q2+w3×Q3+w4×Q4

[0399] Quantum annealing utilizes the quantum tunneling effect to escape local optima and find the global optimum, thus transforming the aforementioned QUBO matrix Q... combined Mapped to the Ising model,

[0400] E(s) = -Σ ij J ij s i s j -Σ i h i s i

[0401] Where s i ∈{-1,+1} is the spin variable, converting s to binary q, q i =(s i +1) / 2

[0402] The meaning of this replacement is:

[0403] When spin up s i When =+1, it corresponds to binary q i =1

[0404] When spin down s i When = -1, it corresponds to binary q. i =0;

[0405] Q combined Convert to Ising parameters J and h.

[0406] Decode q into environmental control parameters X, calculate the performance of this scheme on four objectives [f1*,f2*,f3*,f4*], and thus obtain all the solution sets as Pareto optimal solution sets. The Pareto optimal solution set contains 10 schemes, a complete spectrum from "healthiest but most energy-intensive" to "most energy-efficient but slightly detrimental to health", and users can choose according to their current needs.

[0407] The quantum annealing algorithm described above successfully transforms complex multi-objective optimization problems into forms that can be efficiently solved by quantum computing, and generates diverse Pareto optimal solution sets, providing users with scientific and flexible environmental control solutions.

[0408] In step S160, the environment and equipment clusters are regulated based on the Pareto optimal solution set and the biological rhythm digital twin model.

[0409] In this embodiment, the strategy for selecting the final control scheme from the Pareto optimal solution set obtained through the above steps is as follows:

[0410] The utility function U for each Pareto optimal solution is calculated as follows:

[0411] U = α·C + (1-α)·E

[0412] Where: C is the normalized comfort index, with a value ranging from 0 to 1; E is the normalized energy efficiency index, with a value ranging from 0 to 1; α is the weighting parameter, with a value ranging from 0 to 1.

[0413] The weighting parameter α is dynamically adjusted according to different time periods: During the daytime period (6:00-22:00): α is set to 0.7, prioritizing comfort; During the nighttime period (22:00-6:00): α is set to 0.3, prioritizing energy saving.

[0414] Calculate the utility function values ​​of all Pareto optimal solutions, and select the solution with the largest utility function value as the control scheme for the current moment.

[0415] The Pareto optimal solution set obtained through the above steps covers a complete spectrum of environmental control parameter schemes, ranging from "healthiest but most energy-intensive" to "most energy-efficient but slightly detrimental to health." This can be combined with the physiological state prediction model in the biorhythm digital twin model to predict user needs and select the environmental control parameter scheme that best meets the user's requirements, enabling the user to experience a comfortable and energy-efficient environment.

[0416] In one embodiment of this application, the biorhythm digital twin model includes multi-scale biorhythm decomposition results, a physiological rhythm coupling network, and a physiological state prediction model. The construction of the biorhythm digital twin model involves performing empirical mode decomposition on the physiological state vector.

[0417] Empirical mode decomposition was performed on the physiological state vector to obtain multi-scale biological rhythm decomposition results;

[0418] Granger causality analysis was used to analyze the decomposition results of multi-scale biological rhythms, and the Granger causality analysis results were obtained.

[0419] The transfer entropy results were obtained by calculating the multi-scale biological rhythm decomposition results.

[0420] Based on the results of multi-scale biological rhythm decomposition, the types of influence were determined.

[0421] Based on Granger causality analysis, transfer entropy results, and influence type results, a physiological rhythm coupling network was constructed.

[0422] Based on the multi-scale rhythm decomposition results, the physiological rhythm coupling network, and the physiological state vector, a physiological state prediction model is obtained.

[0423] In the embodiments of this application, physiological state vectors are identified, including heart rate (one data point per minute, for 30 consecutive days, totaling 43,200 data points), heart rate variability (one data point per 5 minutes, for 30 consecutive days, totaling 8,640 data points), respiratory rate (one data point per minute), and body temperature (one data point per 10 minutes). These physiological state data contain complex periodic components and aperiodic noise. Local extreme points in the physiological state data are identified. Taking heart rate variability data as an example, the heart rate variability data collected over the entire time series is scanned to find all local maxima and minima. For example, in a 24-hour segment, 96 maxima and 96 minima may be identified.

[0424] Construct upper and lower envelopes. Using cubic spline interpolation, connect all local maxima to form a smooth upper envelope, and connect all local minima to form a smooth lower envelope. The upper and lower envelopes enclose the heart rate variability data.

[0425] Calculate the envelope mean and extract the first intrinsic mode. Calculate the average of the upper and lower envelopes to obtain the local mean curve. Subtract this mean curve from the heart rate variability data to obtain the first candidate intrinsic mode component. If this component meets the conditions for an intrinsic mode (the difference between the number of zero crossings and the number of extreme points does not exceed 1, and the mean of the envelope is zero everywhere), it is determined as the first intrinsic mode IMF1; otherwise, this component is used as a new input signal, and the above process is repeated until the conditions are met.

[0426] Iterative decomposition yields multiple intrinsic modes (IMFs). The first IMF1 is subtracted from the heart rate variability data to obtain the first residual. This residual is used as the new input signal, and the above three steps are repeated to obtain the second IMF2. This process continues until the residual becomes a monotonic function or has fewer than two extreme points, ultimately yielding the first, second, third, fourth, fifth, sixth, seventh, and eighth IMFs. Spectral analysis is performed on each decomposed IMF to determine its dominant frequency and period. Rhythms are then categorized based on period length.

[0427] Ultra-short rhythms (5-20 minute cycles): usually correspond to IMF1 and IMF2, reflecting rapid regulation of the autonomic nervous system, such as heart rate fluctuations caused by baroreceptor reflexes.

[0428] Short-term rhythms (1-4 hour cycles): usually correspond to IMF3 and IMF4, reflecting sleep cycles, eating rhythms, etc.

[0429] Circadian rhythm (24-hour cycle): usually corresponds to IMF5 or IMF6, reflecting physiological fluctuations controlled by the biological clock.

[0430] Weekly rhythm (7-day cycle): usually corresponds to IMF7 or IMF8, reflecting the difference in lifestyle patterns between weekdays and weekends.

[0431] Multiscale rhythm decomposition results contain a complete feature description of each rhythm component. Taking diurnal rhythms as an example, the output includes: cycle length (e.g., 24.3 hours, indicating an endogenous rhythm slightly longer than 24 hours), amplitude (e.g., the diurnal fluctuation amplitude of heart rate variability is 15 milliseconds), phase (e.g., the peak occurs at 3 pm, delayed by 2 hours relative to the normal population), waveform characteristics (e.g., steep rising edge, gentle falling edge), energy percentage (e.g., accounting for 35% of the total signal energy), and stability index (e.g., a phase standard deviation of 0.5 hours indicates rhythm stability).

[0432] Granger causality analysis is used to determine whether one physiological rhythm has predictive power over another, thereby establishing causal relationships.

[0433] Construct an autoregressive model. For the target rhythm Y (such as a heart rate variability diurnal rhythm), first construct an autoregressive model that uses only its own historical values ​​(heart rate variability diurnal rhythm). The autoregressive model uses values ​​from the past p time points to predict the current value, where the value of p is determined by the Akaike Information Criterion (AIC), typically between 3 and 10. Calculate the prediction error variance of the autoregressive model, denoted as VAR1.

[0434] Construct a joint autoregressive model. Based on the autoregressive model, incorporate historical values ​​of the underlying causal rhythm X (such as the body temperature diurnal rhythm) to construct a joint autoregressive model. This joint autoregressive model uses both historical values ​​of Y and X to predict the current value of Y. Calculate the variance of the prediction error of the joint autoregressive model, denoted as VAR².

[0435] Calculate the Granger causality index. If VAR2 is significantly smaller than VAR1 (significance is determined by the F-test), it indicates that historical information about X helps predict Y, meaning there is a Granger causal relationship between X and Y. The causality strength index is calculated as: GCI = (VAR1 - VAR2) / VAR1, with a value ranging from 0 to 1. The closer the value is to 1, the stronger the causal relationship.

[0436] Determine the causal delay. By testing different time delays (from 1 minute to 60 minutes), identify the delay time that maximizes the Granger causality index. For example, it might be found that changes in body temperature have the strongest effect on heart rate variability after 15 minutes.

[0437] Transfer entropy quantifies the amount of information transferred between rhythms from an information theory perspective, supplementing the linear assumptions of Granger causal analysis.

[0438] State-space discretization. The continuous intrinsic modes in the multi-scale biological rhythm decomposition results are discretized into a finite number of states. For example, the range of heart rate variability diurnal rhythm (30-60 milliseconds) is evenly divided into 10 intervals, each interval representing a state.

[0439] Calculate the conditional entropy. Given the historical states of the target rhythm Y, calculate the uncertainty of the future state of Y (conditional entropy H1). Then, calculate the uncertainty of the future state of Y (conditional entropy H2) given that both the historical states of Y and the historical states of the source rhythm X are known.

[0440] Calculate the transition entropy. Transition entropy TE = H1 - H2, representing the contribution of historical information about X to the prediction of the future state of Y. The unit of transition entropy is bits; a larger value indicates stronger information transfer.

[0441] Determine the type of influence (linear / nonlinear), assume a linear relationship between the two rhythms, and establish a linear regression model. Taking the influence of body temperature diurnal rhythm (X) on heart rate variability diurnal rhythm (Y) as an example, construct a linear model: Y = a × X + b + ε, where a is the slope, b is the intercept, and ε is the error term.

[0442] The parameters a and b are estimated using the least squares method to obtain linear predicted values. Calculate the actual value Y and the predicted value The residual sequence between: If the relationship is indeed linear, the residuals should exhibit a random distribution with a mean of zero and be independent of the independent variable X.

[0443] Residual normality test. The Shapiro-Wilk test is used to determine whether the residuals follow a normal distribution. If the p-value is greater than 0.05, the residuals follow a normal distribution, supporting the linear hypothesis. If the p-value is less than 0.05, the residual distribution deviates from normality, suggesting a possible nonlinear relationship.

[0444] Independence test of residuals. Calculate the autocorrelation function of the residuals to check for correlation between them. If significant autocorrelation exists (autocorrelation coefficient exceeding 0.2), it indicates that the linear model has failed to fully capture the pattern in the data, and there may be nonlinear components.

[0445] Homoscedasticity test of residuals. Divide the range of the independent variable X into 10 intervals and calculate the variance of the residuals in each interval. If the variances of different intervals differ by more than a factor of two, heteroscedasticity is indicated, which is usually a sign of a nonlinear relationship. For example, heart rate variability may fluctuate little at lower body temperatures but fluctuate greatly at higher body temperatures.

[0446] If the residual test of the linear model fails, the system performs more detailed nonlinear feature detection to identify specific nonlinear patterns.

[0447] First, a quadratic term test is performed. A quadratic term is added to the linear model: Y = a1×X + a2×X. 2+b+ε. If the coefficient of the quadratic term a² is significantly non-zero (the p-value of the t-test is less than 0.05), it indicates the existence of a quadratic nonlinear relationship. For example, heart rate variability may exhibit an inverted U-shaped relationship with body temperature, being highest at moderate body temperatures and decreasing at excessively high or low temperatures.

[0448] Next, threshold effect detection is performed. A piecewise linear regression method is used, assuming that the slope of the relationship between rhythms differs before and after a certain threshold point. A grid search method is used to test different potential threshold points (covering the 20% to 80th quantiles of X) to find the threshold that best fits the model. If the goodness of fit of the piecewise model (R²) is high... 2 The improvement of more than 5% compared to a single linear model indicates the existence of a threshold effect. For example, body temperature below 36.5℃ has little effect on heart rate variability, but the effect is significantly enhanced above 36.5℃.

[0449] Third, a saturation effect test is performed. A logarithmic or exponential transformation model is used: Y = a × log(X) + b or Y = a × (1 - exp(-b × X)) + c. If the goodness of fit of the transformed model is significantly improved, it indicates the presence of a saturation effect. For example, a small increase in body temperature can significantly increase heart rate, but when body temperature is already high, further increases tend to saturate the effect on heart rate.

[0450] Mutual information is a correlation measure that does not depend on a specific functional form and can capture any type of dependency, including complex nonlinear relationships.

[0451] First, the two rhythm signals are discretized. The ranges of X and Y are each divided into 20 equally wide intervals. The frequency of each interval combination (Xi,Yj) is counted, and a joint probability distribution P(Xi,Yj) is constructed. At the same time, the marginal probability distributions P(Xi) and P(Yj) are calculated.

[0452] Calculate mutual information value:

[0453]

[0454] If MI(X;Y) is less than the threshold θ (θ takes the value of 0.01), then the edge between nodes X and Y is deleted. Mutual information reflects the degree to which knowing the value of one variable reduces the uncertainty of another variable.

[0455] Compare the relationship between mutual information and the linear correlation coefficient. According to information theory, if two variables are jointly Gaussian distributed (a typical case of linear relationship), there is a definite relationship between mutual information and the correlation coefficient r: MI = -0.5 × log(1 - r). 2 ). Calculate the actual mutual information MI actual And theoretical mutual information MI based on linear correlation coefficient linear If MIactual significantly greater than MI linear (the difference exceeds 0.1 bit), indicating the existence of non - linear components beyond the linear relationship.

[0456] Use rank - based nonparametric methods to further verify the existence of non - linear relationships.

[0457] First, perform Spearman rank - correlation analysis. Convert the original data X and Y into rank data (the positions sorted by size), and calculate the rank - correlation coefficient. If the Spearman correlation coefficient differs from the Pearson linear - correlation coefficient by more than 0.15, it indicates the existence of a monotonic but non - linear relationship. For example, if the Pearson correlation is 0.4 but the Spearman correlation is 0.7, it implies the existence of a monotonically increasing non - linear relationship.

[0458] Secondly, use Kendall's tau test. This method evaluates the consistency of data pairs: the proportion of corresponding Yi < Yj when Xi < Xj. Kendall's tau is sensitive to non - linear monotonic relationships. If Kendall's tau is significant but the R of linear regression 2 is very low (less than 0.3), it strongly implies a non - linear relationship.

[0459] Furthermore, perform distance - correlation analysis. Distance correlation is a method that can detect any type of dependence relationship. Calculate the distance matrix between samples, and then calculate the distance covariance and distance - correlation coefficient. If the distance - correlation coefficient is significantly greater than 0 (the permutation - test p - value is less than 0.05) but the Pearson correlation is close to 0, it indicates the existence of a non - linear relationship.

[0460] Use non - linear machine - learning models to further verify and quantify the strength of non - linear relationships.

[0461] Construct three prediction models: a linear - regression model (baseline), a second - order polynomial - regression model (mildly non - linear), and a random - forest model (capable of capturing any non - linearity). Train these three models using the same training set and evaluate the prediction performance on the test set.

[0462] The performance comparison uses the mean squared error (MSE) of cross - validation. If the MSE of the random forest is more than 20% lower than that of linear regression, it indicates the existence of a significant non - linear relationship. If the MSE of the second - order polynomial is close to that of the random forest (the difference is less than 5%), it indicates that the non - linear relationship is relatively simple and can be approximated by a low - order polynomial. If the random forest is significantly better than the polynomial model, it indicates the existence of complex non - linear patterns.

[0463] Feature importance analysis using random forests is used to identify sources of nonlinearity. High importance of interaction features (different powers of X or products of X with other variables) further confirms the existence of nonlinear relationships.

[0464] Based on the results of the above multiple testing methods, a comprehensive judgment rule is used to determine the final impact type:

[0465] The conditions for determining a linear relationship (all conditions must be met simultaneously):

[0466] R-squared for linear regression 2 Greater than 0.6;

[0467] The residuals passed the normality test (p>0.05);

[0468] The residuals showed no significant autocorrelation (|r|<0.2);

[0469] Homogeneity of residual variance (maximum variance / minimum variance < 2);

[0470] The difference between mutual information and linear expectation is less than 0.1 bits;

[0471] The difference between Pearson correlation and Spearman correlation is less than 0.15;

[0472] The improvement of the random forest model compared to the linear model is less than 10%;

[0473] Cases where a non-linear relationship is determined:

[0474] If any of the above conditions are not met, the relationship is determined to be nonlinear and further classified:

[0475] Quadratic nonlinearity: Quadratic terms are significant and the second-order polynomial model R 2 An increase of more than 10%;

[0476] Threshold nonlinearity: Piecewise regression models significantly outperform linear models and exhibit clear inflection points;

[0477] Saturated nonlinearity: linear fit improved by more than 15% after logarithmic or exponential transformation;

[0478] Complex nonlinearity: Random forests significantly outperform multinomial models and require nonparametric methods for description.

[0479] In coupled networks, each edge is labeled not only with causal strength, information transmission amount, and time delay, but also with detailed indication of its type of influence:

[0480] Example of linear effect annotation: "Linear (slope = 2.3, R 2 =0.75)" indicates a linear relationship of Y = 2.3X + b, with a goodness of fit of 0.75.

[0481] Example of nonlinear effect annotation: "Secondary nonlinearity (dominant term = X)" 2 The inflection point is 36.8℃, indicating a quadratic relationship, with an extreme value appearing near 36.8℃.

[0482] Example of complex nonlinear annotation: "Complex nonlinearity (mutual information = 0.45 bits, table lookup required), indicates complex relationships, requiring querying the stored nonparametric model."

[0483] When constructing a causal model of environmental-physiological relationships, an appropriate modeling method is selected based on the identified type of influence. For linear influences, a simple linear transfer function is used; for nonlinear influences, a corresponding nonlinear function or lookup table method is used. Furthermore, linear relationships can be solved using efficient convex optimization methods; nonlinear relationships require more complex global optimization algorithms, such as quantum annealing or genetic algorithms. Accurate identification of influence types ensures that the entire system uses appropriate mathematical tools, guaranteeing both accuracy and optimized computational efficiency.

[0484] By integrating Granger causality analysis, influence type, and transfer entropy results, a physiological rhythm coupling network is constructed. The network nodes represent different time-scale rhythms of various physiological indicators (e.g., "heart rate-diurnal rhythm," "body temperature-diurnal rhythm," etc.). The directed edges of the network represent the causal relationships between rhythms, and each edge has four attributes: causal strength (Granger causality index), information transfer amount (transfer entropy value), time delay (minutes), and influence type (linear / nonlinear).

[0485] Topological analysis of coupled physiological rhythm networks also identifies key nodes (hubs) and critical paths. For example, if a node with 5 outgoing edges is found to affect the other 5 rhythms, it is marked as a hub node. If a path of "light → melatonin → body temperature → heart rate" is found to form a major information transmission path, it is marked as a critical path.

[0486] In the embodiments of this application, performing Granger causality analysis directly on the physiological state vector will encounter serious confounding problems. For example, the raw heart rate signal contains multiple rhythmic components, and the raw body temperature signal also contains multiple rhythmic components. When analyzing the causal effect of body temperature on heart rate, the following confounding may occur:

[0487] The circadian rhythm (24-hour cycle) of body temperature may affect the circadian rhythm of heart rate.

[0488] Short-term fluctuations in body temperature (2-hour cycle) may affect short-term fluctuations in heart rate.

[0489] Both likely have circadian rhythms controlled by biological clocks, suggesting a common cause rather than a causal relationship.

[0490] The effects of different time scales may mask each other, leading to inaccurate analytical results.

[0491] Therefore, by using the results of multi-scale biological rhythm decomposition to separate rhythms at different time scales, and then conducting causal analysis at the same time scale, the aforementioned confusion problem can be avoided, enabling a more accurate determination of the actual causal relationship of users and effective regulation tailored to users.

[0492] The object of physiological rhythm coupling network analysis is not the physiological state vector, but rather the individual rhythmic components of the multi-scale biological rhythm decomposition result. This is because rhythms at different time scales may have different coupling relationships. For example, body temperature and heart rate may be strongly coupled at the diurnal rhythm scale (both regulated by the biological clock), but relatively independent at the ultrashort rhythm scale (controlled by different rapid regulatory mechanisms). By performing coupling analysis separately for each time scale, a more accurate map of inter-rhythmic relationships can be obtained.

[0493] In the embodiments of this application, the construction process of the Long Short-Term Memory (LSTM) neural network is as follows:

[0494] The results of multi-scale rhythm decomposition are transformed into input features for a neural network. The input vector at each time point includes: the current physiological state value (e.g., heart rate 72 bpm), the current phase of each rhythm component (e.g., diurnal rhythm phase π / 4), the expected contribution of each rhythm component (e.g., diurnal rhythm contribution +5 bpm), and a predicted coupling effect (e.g., a decrease in body temperature is expected to lead to a decrease in heart rate of 3 bpm). The dimension of the input vector is typically between 50 and 100.

[0495] The system comprises an input layer, a feature extraction layer, an attention layer, a fusion layer, and a decoding layer. The input layer includes channels for raw values, rhythm decomposition, feature parameters, coupling effects, and temporal encoding. The feature extraction layer includes ultra-short rhythm layers, short-term rhythm layers, diurnal rhythm layers, and long-term rhythm layers. The attention layer includes an immediate attention head, a short-term attention head, a medium-term attention head, and a long-term attention head. The decoding layer includes a first fully connected layer and a second fully connected layer. Each of the ultra-short rhythm layer, short-term rhythm layer, diurnal rhythm layer, and long-term rhythm layer contains 128 LSTM units. Each LSTM unit includes three gating mechanisms: a forget gate, an input gate, and an output gate. The forget gate determines how much historical information is retained, the input gate determines how much new information is received, and the output gate determines what information is output. This gating mechanism enables LSTM to learn long-term dependencies, making it particularly suitable for processing physiological signals with multi-timescale characteristics.

[0496] Immediate attention heads, short-term attention heads, medium-term attention heads, and long-term attention heads enable the model to automatically focus on the historical time points and features most important for the prediction. Attention weights are learned and reflect the importance of different historical information to the current prediction. For example, when predicting nighttime sleep status, the model automatically increases the attention weight of data from the same period of the previous night.

[0497] The 30-day historical data was divided into a training set (first 25 days), a validation set (days 26-28), and a test set (days 29-30). Each sample included physiological state data from the past 24 hours as input and physiological state data from the next 8 hours as output labels. A combined loss function was used, including: mean squared error loss (for continuous value prediction, such as heart rate and body temperature), cross-entropy loss (for classification prediction, such as sleep stage), and phase consistency loss (to ensure the predictor's rhythm phase is reasonable). The total loss was a weighted sum of the three, with weights set according to the importance of the prediction task.

[0498] The Adam optimizer is used with an initial learning rate of 0.001 and a batch size of 32. During training, the validation set loss is monitored; the learning rate is halved if the validation set loss does not decrease for 5 consecutive epochs. Training is stopped if the validation set loss does not decrease for 10 consecutive epochs to prevent overfitting. Training typically takes 50-100 epochs. Multiple LSTM models with different initialization parameters (usually 5) are trained, and their predictions are integrated using a weighted average to improve the stability and accuracy of the predictions. The weights are determined based on the performance of each model on the validation set, resulting in a physiological state prediction model. This model can receive the current physiological state and historical 24-hour physiological data, and output physiological state predictions for the next 8 hours (32 time points, one every 15 minutes). The predictions for each time point include: heart rate (beats / min), heart rate variability SDNN (milliseconds), heart rate variability RMSSD (milliseconds), core body temperature (degrees Celsius), sleep stage (awake / N1 / N2 / N3 / REM), sleep stage probability distribution, alertness level (0-1 continuous value), and prediction confidence (0-1, based on model uncertainty estimation).

[0499] The physiological state prediction model also identifies key physiological events, including: expected sleep onset time (the time when the sleep stage transitions from wakefulness to N1), deep sleep window (the period when the N3 stage lasts for more than 20 minutes), REM sleep cycle (the start and end times of the REM stage), and wake-up time (the time of the last transition from sleep to wakefulness). Each event is accompanied by a probability of occurrence, which helps subsequent steps assess the reliability of the prediction.

[0500] In one embodiment of this application, empirical mode decomposition is performed on the physiological state vector to obtain multi-scale biological rhythm decomposition results, including:

[0501] Identify the physiological state vector to obtain all local maxima and all local minima;

[0502] Using cubic spline interpolation, all local maxima are connected to form a smooth upper envelope and all local minima are connected to form a smooth lower envelope. The upper and lower envelopes enclose the physiological state vector.

[0503] Calculate the average value of the upper and lower envelopes to obtain the local mean curve;

[0504] Subtracting the local mean curve from the physiological state vector yields the first candidate intrinsic mode component.

[0505] Determine whether the first candidate intrinsic mode component satisfies the intrinsic mode condition;

[0506] If the first candidate intrinsic mode component satisfies the intrinsic mode condition, it is determined to be the first intrinsic mode.

[0507] If the first candidate intrinsic mode component does not meet the intrinsic mode condition, the first candidate intrinsic mode component is used as the physiological state vector, and the above process is repeated until the first candidate intrinsic mode component meets the intrinsic mode condition, thus obtaining the first intrinsic mode.

[0508] The first residual is obtained by subtracting the first eigenmode from the physiological state vector;

[0509] Using the first residual as the physiological state vector, repeat the above process until the obtained first residual becomes a monotonic function or has fewer than 2 extreme points. Then, the first eigenmode in the repeated process is successively determined as the second eigenmode, the third eigenmode, the fourth eigenmode, the fifth eigenmode, the sixth eigenmode, the seventh eigenmode, and the eighth eigenmode.

[0510] Based on the first, second, third, fourth, fifth, sixth, seventh, and eighth intrinsic modes, the multi-scale biological rhythm decomposition results are obtained.

[0511] In the embodiments of this application, physiological state vectors are identified, including heart rate (one data point per minute, for 30 consecutive days, totaling 43,200 data points), heart rate variability (one data point per 5 minutes, for 30 consecutive days, totaling 8,640 data points), respiratory rate (one data point per minute), and body temperature (one data point per 10 minutes). These physiological state data contain complex periodic components and aperiodic noise. Local extrema of the physiological state data are identified. Taking heart rate variability data as an example, the heart rate variability data collected over the entire time series is scanned to find all local maxima and minima. For example, in a 24-hour segment, 96 maxima and 96 minima may be identified. Assume that in the first hour (12 data points), 5 local maxima and 5 local minima are found. The locations and values ​​of these extrema are as follows:

[0512] Maximum points: 2nd point (42ms), 5th point (45ms), 7th point (43ms), 9th point (44ms), 11th point (46ms)

[0513] Minimum points: 1st point (38ms), 4th point (37ms), 6th point (39ms), 8th point (38ms), 10th point (36ms)

[0514] Construct upper and lower envelopes. Using cubic spline interpolation, connect all local maxima to form a smooth upper envelope, and connect all local minima to form a smooth lower envelope. The upper and lower envelopes enclose the heart rate variability data.

[0515] Calculate the envelope mean and extract the first intrinsic mode. Calculate the average of the upper and lower envelopes to obtain the local mean curve. Subtract this mean curve from the heart rate variability data to obtain the first candidate intrinsic mode component. If this component meets the conditions for an intrinsic mode (the difference between the number of zero crossings and the number of extreme points does not exceed 1, and the mean of the envelope is zero everywhere), it is determined as the first intrinsic mode IMF1; otherwise, this component is used as a new input signal, and the above process is repeated until the conditions are met.

[0516] Iterative decomposition yields multiple intrinsic modes. The first intrinsic mode (IMF1) is subtracted from the heart rate variability data to obtain the first residual, r1. This residual is used as the new input signal, and the processes of extreme point identification, envelope construction, and mean calculation are repeated. Since the highest frequency components have been removed, the number of extreme points in r1 is significantly reduced. Within the same hour, there may be only 3 maxima and 3 minima, corresponding to a period of approximately 20-30 minutes. After screening, IMF2 is obtained, representing a slightly slower physiological rhythm, such as the Mayer wave in blood pressure. This process continues until the residual becomes a monotonic function or has fewer than 2 extreme points, ultimately yielding the first, second, third, fourth, fifth, sixth, seventh, and eighth intrinsic modes. Spectral analysis is performed on each decomposed intrinsic mode to determine its dominant frequency and period. Rhythms are categorized based on period length.

[0517] Ultra-short rhythms (5-20 minute cycles): usually correspond to IMF1 and IMF2, reflecting rapid regulation of the autonomic nervous system, such as heart rate fluctuations caused by baroreceptor reflexes.

[0518] Short-term rhythms (1-4 hour cycles): usually correspond to IMF3 and IMF4, reflecting sleep cycles, eating rhythms, etc.

[0519] Circadian rhythm (24-hour cycle): usually corresponds to IMF5 or IMF6, reflecting physiological fluctuations controlled by the biological clock.

[0520] Weekly rhythm (7-day cycle): usually corresponds to IMF7 or IMF8, reflecting the difference in lifestyle patterns between weekdays and weekends.

[0521] Multiscale rhythm decomposition (EMD) results contain a complete feature description of each rhythm component. Taking diurnal rhythms as an example, the output includes: cycle length (e.g., 24.3 hours, indicating an endogenous rhythm slightly longer than 24 hours), amplitude (e.g., the diurnal fluctuation amplitude of heart rate variability is 15 milliseconds), phase (e.g., the peak occurs at 3 PM, delayed by 2 hours relative to the normal population), waveform characteristics (e.g., steep rising edge, gentle falling edge), energy percentage (e.g., accounting for 35% of the total signal energy), and stability index (e.g., a phase standard deviation of 0.5 hours, indicating rhythm stability). A key feature of the EMD method is its adaptive decomposition based on the inherent characteristics of the signal, rather than a preset frequency. This makes it particularly suitable for analyzing non-stationary physiological signals. Each IMF is arranged naturally from high frequency to low frequency: IMF1 and IMF2 contain the highest frequency components (5-20 minute cycles), which correspond precisely to the rapid regulation of the autonomic nervous system; IMF3 and IMF4 contain the medium frequency components (1-4 hour cycles), which correspond to sleep stage transitions; IMF5 and IMF6 contain the low frequency components (12-24 hour cycles), which correspond to circadian rhythms; and IMF7 and IMF8 contain the lowest frequency components (3-7 day cycles), which correspond to the periodic changes in lifestyle habits.

[0522] In this process, a Hilbert Transform is performed on each IMF to obtain the instantaneous frequency and instantaneous phase.

[0523] Amplitude reflects the intensity of rhythm fluctuations and is obtained through envelope analysis. Identify all local maxima in the IMF. For IMF6 (circadian rhythm), there are approximately 30 major maxima over 30 days (one per day). Calculate the envelope amplitude. For each cycle, amplitude = (local maxima - local minima) / 2. Calculate the mean amplitude and amplitude variability. Mean amplitude: the average amplitude across all cycles. For example, the mean amplitude of a heart rate variability circadian rhythm is 15 ms, meaning the variation from lowest to highest amplitude is approximately 30 ms. Amplitude standard deviation: reflects the stability of the amplitude. A smaller standard deviation indicates a stable rhythm intensity. Calculate the temporal trend of amplitude variation. Use linear regression analysis to analyze amplitude changes over 30 days to determine whether the rhythm is strengthening or weakening.

[0524] Phase tells us when the rhythm reaches its peak or trough during the day. Determine a reference time point. Midnight (0:00) is usually chosen as the reference point for phase 0. Identify peak times. For IMF6, find the time of the daily maximum. For example, the peak times for 30 consecutive days are: 15:10, 15:15, 15:08, 15:20...; Calculate the average phase. Convert time to phase angle: the phase corresponding to 15:00 = (15 / 24) × 2π = 5π / 4 radians. Calculate the 30-day average phase. Compare with the standard phase. Heart rate variability in the normal population is usually highest between 4-6 AM (corresponding to the peak of parasympathetic activity). If an individual's heart rate peaks at 3 PM, it indicates a phase delay of approximately 9 hours, which may suggest a circadian rhythm disorder. Calculate phase stability. The standard deviation of phase = 0.5 hours represents the diurnal variability of peak time. The smaller the standard deviation, the more stable the rhythm.

[0525] Waveform characteristics describe the shape of the rhythm, unlike a simple sine wave. Calculate rise and fall times. Rise time: the time required from the trough to the peak; Fall time: the time required from the peak to the trough. Calculate waveform skewness. Skewness measures the asymmetry of the waveform: Skewness > 0: the waveform is right-skewed, with a steep rise and a gentle fall; Skewness < 0: the waveform is left-skewed, with a gentle rise and a steep fall; Skewness = 0: a symmetrical waveform, similar to a sine wave. The third step is to calculate waveform kurtosis. Kurtosis describes the sharpness of the peaks: Kurtosis > 3: the peaks are sharp, with changes concentrated in a short time; Kurtosis < 3: the peaks are flat, with changes dispersed over a longer time. Harmonic analysis. Fourier analysis determines the proportion of each harmonic in the waveform. A pure sine wave only has a fundamental frequency, while actual physiological rhythms usually contain multiple harmonics, causing the waveform to deviate from a sine wave.

[0526] Energy percentage reflects the contribution of each rhythmic component to the overall signal variation. Calculate the energy of each IMF. For discrete signals, energy E = Σ(IMF value). 2 .For example:

[0527] IMF1 energy: E1 = 12000 (unit: ms) 2 IMF2 energy: E2 = 8000; IMF3 energy: E3 = 5000; IMF4 energy: E4 = 6000; IMF5 energy: E5 = 4000; IMF6 energy: E6 = 15000 (circadian rhythm energy is the largest); IMF7 energy: E7 = 3000; IMF8 energy: E8 = 2000; residual energy: E9 = 1000.

[0528] Calculate the total energy: E total =E1+E2+...+E9=56000

[0529] Calculate the energy percentage of each component. IMF6 (circadian rhythm) energy percentage = 15000 / 56000 = 26.8%; this means that the circadian rhythm explains 26.8% of the total variation in heart rate variability.

[0530] Explanation of the physiological significance of energy percentage: High energy percentage of circadian rhythm (>30%): indicates strong biological clock regulation and obvious circadian rhythm; High energy percentage of ultrashort rhythm (>40%): may indicate autonomic nervous system dysfunction; High energy percentage of residual (>20%): the signal contains more random components or unidentified rhythms.

[0531] Stability reflects the temporal consistency of rhythms and includes multiple dimensions:

[0532] Phase stability. Calculate the standard deviation of the peak time each day over 30 days. A standard deviation of 0.5 hours indicates that the peak time is very stable, with diurnal variation of less than 30 minutes.

[0533] Amplitude stability. Calculate the coefficient of variation (standard deviation / mean) of the amplitude. A coefficient of variation <0.2 indicates stable amplitude; >0.5 indicates large amplitude variations and unstable rhythm.

[0534] Periodic stability. Calculate the standard deviation of the instantaneous period. For diurnal rhythms, if the standard deviation of the period is <0.5 hours, the rhythm period is considered very stable.

[0535] Phase coherence. Calculate the phase consistency between adjacent cycles. Quantize using Phase Lock Value (PLV): PLV close to 1 indicates high phase consistency, and close to 0 indicates random phase variation.

[0536] Overall stability index. This integrates the above four indicators:

[0537] Stability index = 0.3 × phase stability + 0.3 × amplitude stability + 0.2 × period stability + 0.2 × phase coherence;

[0538] Stability index > 0.8: The rhythm is very stable and the physiological regulatory function is good;

[0539] Stability index 0.5-0.8: The rhythm is basically stable, but may be subject to slight disturbances;

[0540] Stability index <0.5: The rhythm is unstable and there may be physiological dysfunction.

[0541] In one embodiment of this application, Granger causality analysis is used to analyze the decomposition results of multi-scale biological rhythms, and the Granger causality analysis results include:

[0542] Granger causality analysis was used to pair the results of multi-scale biological rhythm decomposition based on the principle of similar cycles to obtain the causal relationships between various biological rhythms.

[0543] Based on the multi-scale biological rhythm decomposition results, a first autoregressive model was constructed.

[0544] By incorporating the causal relationships between corresponding biological rhythms into the first autoregressive model, the second autoregressive model is determined.

[0545] Determine the causal strength based on the first and second autoregressive models;

[0546] Determine the delay time based on the strength of causality;

[0547] Based on the causal strength and the delay time, the Granger causality analysis results are obtained.

[0548] In this embodiment, the multi-scale circadian rhythm decomposition results have yielded multiple intrinsic mode functions (IMFs) for each physiological indicator. Now, it is necessary to pair the IMFs of the same timescale for different physiological indicators. Specific pairing method:

[0549] Taking heart rate and body temperature as examples:

[0550] Heart rate IMF1 (cycle 12 minutes), IMF2 (cycle 25 minutes), ..., IMF6 (cycle 24.6 hours)

[0551] IMF1 (15-minute cycle), IMF2 (28-minute cycle), ..., IMF6 (24.2-hour cycle) for body temperature.

[0552] The matching principle is to pair IMFs with similar operating cycles:

[0553] Ultra-short rhythm pairing: heart rate first intrinsic modality IMF1 (12 minutes) and body temperature first intrinsic modality IMF1 (15 minutes);

[0554] Short-term rhythm pairing: heart rate third intrinsic modality IMF3 (90 minutes) and body temperature third intrinsic modality IMF3 (85 minutes);

[0555] Circadian rhythm pairing: heart rate sixth intrinsic modality IMF6 (24.6 hours) and body temperature sixth intrinsic modality IMF6 (24.2 hours);

[0556] Each IMF pair undergoes an independent Granger causality analysis to determine the causal relationship at that time scale.

[0557] Construct an autoregressive model for each pair of IMFs

[0558] Taking the analysis of circadian rhythms (heart rate IMF6 and body temperature IMF6) as an example, the specific process is as follows:

[0559] The first autoregressive model of heart rate circadian rhythm uses historical values ​​of heart rate IMF6 to predict current values. IMF6 is a rhythmic component with a 24-hour cycle, with a representative value taken every hour, totaling 720 data points over 30 days.

[0560] Model format: HR IMF6 (t)=a1×HR IMF6 (t-1)+a2×HR IMF6 (t-2)+...+a p ×HR IMF6 (tp)+ε t (t)

[0561] in:

[0562] HR IMF6 (t) is the value of the diurnal rhythm of heart rate at time t; HR IMF6 (t-1) is the value from 1 hour ago, HR IMF6 (t-2) is the value from 2 hours ago, p is the model order, determined by the Akaike Information Criterion (AIC), assuming p = 6 (using data from the past 6 hours); a1 to a p ε1(t) is the regression coefficient, estimated using the least squares method; ε1(t) is the prediction error.

[0563] Calculation process:

[0564] Using the first 25 days (600 data points) as training data, the following values ​​were obtained using the least squares method: a1 = 0.8, a2 = 0.3, a3 = -0.1, a4 = -0.05, a5 = 0.02, a6 = 0.01;

[0565] The prediction error variance VAR1 is calculated to be 2.5 (unit: predictions). 2 / point 2 ).

[0566] Including historical IMF6 body temperature values ​​in the second autoregressive model to see if it improves predictions; the second autoregressive model takes the following form:

[0567] HR IMF6 (t)=a1×HR IMF6 (t-1)+...+a6×HR IMF6 (t-6)+b1×T -IMF6 (t-1)+b2×T -IMF6 (t-2)+…+b6×T -IMF6 (t-6)+ε2(t)

[0568] Wherein: TIMF6 (tk) is the value of the body temperature diurnal rhythm k hours ago; b1 to b6 are the coefficients of the influence of body temperature on heart rate.

[0569] Calculation process: Using the same training data, new coefficients are obtained through the least squares method, with slight adjustments to the original coefficients: a1 = 0.75, a2 = 0.28...

[0570] Correlation coefficients of body temperature: b1 = 0.15, b2 = 0.25, b3 = 0.20, b4 = 0.10, b5 = 0.05, b6 = 0.02

[0571] Note that b2 is the largest, indicating that the body temperature 2 hours ago has the greatest impact on the current heart rate; calculate the new prediction error variance VAR2 = 1.8;

[0572] The causal strength is calculated based on the prediction error variance of the two models:

[0573] Granger causality index GCI = (VAR1 - VAR2) / VAR1 = (2.5 - 1.8) / 2.5 = 0.28;

[0574] This means that after incorporating information on the body temperature circadian rhythm, the prediction error of the heart rate circadian rhythm was reduced by 28%.

[0575] Perform a statistical significance test (F-test):

[0576] F-statistic = ((VAR1-VAR2) / 6) / (VAR2 / (600-12)) = (0.7 / 6) / (1.8 / 588) = 38.2;

[0577] At a significance level of 0.05, the critical value of F is approximately 2.1;

[0578] Since 38.2 > 2.1, the effect of body temperature diurnal rhythm on heart rate diurnal rhythm is statistically significant.

[0579] By testing different delay times, the strongest causal relationship can be identified:

[0580] Tests were conducted for delays ranging from 1 to 6 hours:

[0581] Delay by 1 hour: Use only T_IMF6(t-1), GCI = 0.08

[0582] Delay by 2 hours: Use only T_IMF6(t-2), GCI = 0.18

[0583] Delay by 3 hours: Use only T_IMF6(t-3), GCI = 0.12

[0584] Delay by 4 hours: Use only T_IMF6(t-4), GCI = 0.06

[0585] The results showed that the circadian rhythm of body temperature had the strongest effect on the circadian rhythm of heart rate after 2 hours. This is consistent with physiological knowledge: changes in body temperature affect heart rate by influencing metabolic rate, a process that takes time.

[0586] For example, Granger causality analysis was performed on heart rate IMF1 (12-minute cycle) and body temperature IMF1 (15-minute cycle): Due to the very short cycle, minute-level data was used, with a model order p=10 (using the past 10 minutes). The VAR1 of the heart rate model alone was 0.8; the VAR2 after adding body temperature was 0.75; the GCI was 0.063 (weak causality); the optimal delay was 1-2 minutes. On ultra-short timescales, the direct effect of body temperature on heart rate is very small (GCI only 0.063). This is reasonable because this timescale mainly reflects rapid physiological regulation such as respiration, and changes in body temperature are too slow to have a significant effect.

[0587] Causal analysis at the short-term circadian rhythm level (IMF3 and IMF4) analyzed heart rate IMF3 (90-minute cycle) and body temperature IMF3 (85-minute cycle): using data at 15-minute intervals, model order p=8 (using the past 2 hours): VAR1 for the heart rate model alone was 1.5, VAR2 with body temperature included was 1.2, GCI was 0.20 (moderate causality), and the optimal delay was 15-30 minutes. On the sleep cycle timescale, body temperature rhythm has a moderate effect on heart rate rhythm. This reflects the physiological phenomenon that a decrease in body temperature during sleep leads to a decrease in heart rate.

[0588] Causal analysis at the circadian rhythm level (IMF5 and IMF6) is the most important level of analysis, as described in detail above: GCI = 0.28 (strong causal relationship), optimal delay: 2 hours; this reflects the regulatory effect of circadian body temperature rhythm on heart rate rhythm.

[0589] Causal analysis at the weekly rhythm level (IMF7 and IMF8) was performed on heart rate IMF8 (7.2-day cycle) and body temperature IMF8 (7.0-day cycle): using daily average data, model order p=7 (using the past week): VAR1 of the heart rate model alone was 0.3, VAR2 after adding body temperature was 0.28, GCI was 0.067 (very weak causality). The weekly rhythm mainly reflects lifestyle habits (weekdays vs. weekends). The weekly rhythms of body temperature and heart rate may both be the result of lifestyle, rather than mutually causal.

[0590] Constructing a multi-scale causal relationship map: Integrating Granger causal analysis results from all time scales to form a complete description of causal relationships:

[0591] Multiscale effects of body temperature on heart rate: Ultra-short rhythm (10-20 minutes): GCI = 0.063, delayed by 1-2 minutes, very weak effect; Short-duration rhythm (1-4 hours): GCI = 0.20, delayed by 15-30 minutes, moderate effect; Circadian rhythm (24 hours): GCI = 0.28, delayed by 2 hours, strong effect; Weekly rhythm (7 days): GCI = 0.067, delayed by 1 day, very weak effect. This multiscale analysis reveals an important finding: body temperature mainly affects heart rate through circadian and short-duration rhythms, while having a very small effect on the ultra-short-weekly rhythm scale.

[0592] Analysis using rhythm decomposition clearly identified strong causal relationships at the diurnal rhythm level (GCI = 0.28), accurately located the optimal delay of 2 hours, and found that there were almost no causal relationships at the ultrashort rhythm level, providing precise time-scale information for subsequent prediction and control.

[0593] In one embodiment of this application, a physiological state prediction model is obtained based on the multi-scale rhythm decomposition results, the physiological rhythm coupling network, and the physiological state vector, including:

[0594] Based on the results of multi-scale rhythm decomposition and physiological rhythm coupling network, a first long short-term memory network model is constructed. The network model includes an input layer, a feature extraction layer, an attention layer, a fusion layer, and a decoding layer. The input layer includes a raw value channel, a rhythm decomposition channel, a feature parameter channel, a coupling effect channel, and a time encoding channel. The feature extraction layer includes ultra-short rhythm layering, short-term rhythm layering, diurnal rhythm layering, and long-term rhythm layering. The attention layer includes an immediate attention head, a short-term attention head, a medium-term attention head, and a long-term attention head. The decoding layer includes a first fully connected layer and a second fully connected layer.

[0595] The first long short-term memory network model is initialized to obtain the initial long short-term memory network model;

[0596] The initial long short-term memory network model is trained and tested based on the physiological state vector to obtain the physiological state prediction model.

[0597] In the embodiments of this application, the input layer of the neural network is not a simple channel for direct input of physiological state data, but a multi-channel data structure constructed based on the multi-scale rhythm decomposition results and the physiological rhythm coupling network. Each channel carries a specific type of information, ensuring that the network can fully understand various aspects of the physiological state. Specifically:

[0598] The first input channel is the raw value channel, which takes in the user's physiological data, including heart rate, body temperature, systolic blood pressure, diastolic blood pressure, and heart rate variability. For example, at 9 PM, the input vector would be [75 beats per minute, 36.8 degrees Celsius, 125 mmHg, 82 mmHg, 35 milliseconds]. These raw values ​​provide a baseline level of physiological function. However, raw values ​​alone are insufficient because a heart rate of 75 beats per minute could be in a rising phase (morning) or a falling phase (evening), which the network cannot distinguish. Therefore, a second channel is needed.

[0599] The second input channel is the rhythm decomposition channel, which takes into account the current values ​​of each IMF component obtained from the multi-scale rhythm decomposition. For heart rate, this includes the values ​​of eight IMFs, for example, [2.3, -1.5, 0.8, -0.3, 1.2, 4.5, 0.2, -0.1]. Each value represents the current contribution of a specific periodic rhythm. The first value, 2.3, is the current value of IMF1, representing a current contribution of 2.3 heart rate increases per minute from the ultrashort rhythm with a 12-minute cycle; the sixth value, 4.5, is the current value of IMF6, representing a current contribution of 4.5 heart rate increases per minute from the diurnal rhythm with a 24.6-hour cycle. Through these decomposed values, the network knows that the current heart rate of 75 is composed of the baseline value of 68 plus the contributions of each rhythm.

[0600] The third input channel is the feature parameter channel. Knowing only the IMF value is insufficient; the network also needs to understand the dynamic characteristics of each rhythm. This channel provides four key features for each IMF: instantaneous period, instantaneous phase, instantaneous amplitude, and phase stability. Taking the circadian rhythm of heart rate (IMF6) as an example, the input feature vector is [24.6 hours, 7π / 4 radians, 8.5 beats per minute, 0.92]. This tells the network that the rhythm's period is 24.6 hours (slightly longer than 24 hours), the current phase is 7π / 4 (equivalent to around 10 PM), the amplitude is 8.5 beats per minute (the range of change from lowest to highest), and the phase stability of 0.92 indicates that the rhythm is very regular.

[0601] The fourth input channel is the coupling effect channel, which takes into account the interactions between physiological indicators calculated by the circadian rhythm coupling network. For predicting heart rate, the expected impact of other indicators on heart rate needs to be input. For example, the influence vector of body temperature on heart rate is [0.2, 0.8, 0.5, 0.3, 0.1, 0.05], where these six values ​​represent the intensity of the influence of body temperature changes on heart rate over the next 1 to 6 hours. The second value, 0.8, is the largest, indicating that body temperature changes mainly affect heart rate after 2 hours. If the current body temperature drops by 0.5 degrees Celsius, the network can calculate that the heart rate will decrease by 0.5 × 3.2 × 0.28 = 0.45 beats per minute after 2 hours, where 3.2 is the body temperature-heart rate influence coefficient and 0.28 is the Granger causality index.

[0602] The fifth input channel is the time-encoded channel. The circadian rhythm is closely related to time; the network needs to know "what time it is now." This channel's input includes: current hour (0-23), day of the week (1-7), date (1-31), month (1-12), season code (1-4), whether it's a weekday (0 or 1), and whether it's a holiday (0 or 1). For example, [21,3,15,10,3,1,0] represents 9 PM, Wednesday, October 15th, autumn, weekday, and non-holiday. This information helps the network understand the current time context and is particularly important for predicting weekly rhythms.

[0603] Standard LSTM processes time series data through three gating mechanisms (input gate, forget gate, and output gate) and a memory unit. The input gate determines which information should be remembered, the forget gate determines which previous memories should be forgotten, the output gate determines what information should be output, and the memory unit stores long-term information. However, standard LSTM treats information at all time scales equally, which is unsuitable for predicting physiological signals. We need to inform the network that ultrashort rhythms (minute-level) should be updated and forgotten quickly, while diurnal rhythms (24-hour) should remain stable and updated slowly, and the prediction methods for different rhythms should also differ. Therefore, the feature extraction layer (LSTM layer) in this embodiment includes ultrashort rhythm layers, short-term rhythm layers, diurnal rhythm layers, and long-term rhythm layers.

[0604] The first layer processes ultrashort rhythm hierarchies. The input is sequences from IMF1 and IMF2, with each time step having 10 dimensions (5 indicators × 2 IMFs). The sequence length is 60 (60 minutes of data). Including the 4-dimensional time encoding, the total input dimension is 14.

[0605] Create LSTM units: Set the hidden state dimension to 32, which is sufficient to capture ultra-short rhythmic patterns while avoiding overfitting. The memory unit is also 32-dimensional.

[0606] Define the weight matrix:

[0607] Input gate weights Wi: Size (14+32)×32=46×32=1472 parameters;

[0608] Forget gate weight Wf: 46 × 32 = 1472 parameters;

[0609] Candidate memory weights Wc: 46 × 32 = 1472 parameters;

[0610] Output gate weights Wo: 46 × 32 = 1472 parameters;

[0611] Bias vector: 4 × 32 = 128 parameters;

[0612] Total number of parameters: 6016.

[0613] The first-level ultrashort rhythm hierarchical calculation process includes:

[0614] For each time step t (from 1 to 60), the current IMF value (10-dimensional) and the time code (4-dimensional) are concatenated to obtain a 14-dimensional input vector xt. The previous hidden state ht-1 (32-dimensional) and the current input xt (14-dimensional) are concatenated to form a 46-dimensional vector. Multiplied by the input gate weight matrix Wi and with the bias bi, a 32-dimensional vector is obtained. The input gate value it is mapped to the range 0-1 using the sigmoid function (1 / (1+exp(-x))). Similar to the input gate, but with an additional decay coefficient of 0.3, because ultrashort rhythms need to be forgotten quickly. ft = 0.3×sigmoid(Wf×[ht-1,xt]+bf). The tanh activation function (range -1 to 1) is used:

[0615]

[0616] Combining forgetting and input: Because of the small forgetting threshold, old memories decay rapidly.

[0617] Calculate the output gate and hidden state. ot = sigmoid(Wo × [ht-1, xt] + bo), ht = ot × tanh(Ct).

[0618] After 60 time steps, the final hidden state h60 contains the characteristic representation of the ultrashort rhythm.

[0619] The second layer of short-term rhythm hierarchy is constructed. This second layer processes short-term rhythms (IMF3-4) with 64 hidden dimensions and a sequence length of 24 (6 hours, 15-minute intervals). The input includes not only IMF3-4 but also the output of the first layer. Specifically: IMF3-4 sequence (10 dimensions) + first-layer output (32 dimensions) + time encoding (4 dimensions) = 46-dimensional input.

[0620] The weight matrix increases accordingly:

[0621] The weight of each door: (46+64)×64=110×64=7040 parameters;

[0622] Four gates plus offset: 4×7040 + 4×64 = 28416 parameters;

[0623] The forgetting gate coefficient is set to 0.6 to maintain moderate memory. The calculation process is similar to the first layer, but the time step is changed to 15 minutes, resulting in greater variation between adjacent time steps.

[0624] The third layer, the circadian rhythm layer, is constructed. This third layer (IMF5-6) is the most important layer, with 128 hidden dimensions. The input integrates the information from the first two layers: IMF5-6 (10 dimensions) + first layer output (32 dimensions) + second layer output (64 dimensions) + time encoding (4 dimensions) = 110 dimensions.

[0625] With a sequence length of 48 (48 hours of hourly data), it can capture the complete circadian rhythm cycle.

[0626] Weighting scale: (110+128)×128=238×128=30464 parameters per gate, totaling approximately 120,000 parameters. The forgetting gate coefficient is 0.85, preserving long-term memory. This ensures the stable maintenance of the circadian rhythm pattern.

[0627] The fourth layer, a long-term rhythm hierarchy, is constructed to process long-term rhythm layers (IMF7-8 and residuals), with 64 hidden dimensions. The input includes all previous layer outputs: IMF7-8 + residuals (15 dimensions) + the outputs of the first three layers (32 + 64 + 128 = 224 dimensions) + time encoding (4 dimensions) = 243 dimensions. The sequence length is 14 days (14 days of daily data), capturing weekly rhythms. The forgetting gate coefficient is 0.95, almost completely preserving historical information. This layer primarily learns daily life patterns, such as the differences between weekdays and weekends.

[0628] The rhythm-aware capability is embedded into the core computation of the LSTM feature extraction layer. The standard LSTM input gate formula is: it = sigmoid(Wi × [ht-1, xt] + bi).

[0629] After adjusting the LSTM feature extraction layer to rhythm-aware, the calculation is performed in five parts:

[0630] The first part is the contribution of the raw value: Wi_raw × the original measurement value vector; the second part is the contribution of the rhythm state: Wi_rhythm × the current value vector of the IMF; the third part is the contribution of the phase sensing: Wi_phase × the phase encoding vector; the fourth part is the contribution of the coupling effect: Wi_couple × the coupling effect vector; the fifth part is the contribution of the historical information: Wi_hist × the hidden state at the previous time step; the final input gate is: it = sigmoid(Σ five parts + bi) × phase modulation coefficient.

[0631] The phase modulation coefficient is dynamically adjusted according to the current phase. For example, during the rising phase of the circadian rhythm (morning), the coefficient is larger (e.g., 1.2) to encourage the reception of new information; during the steady phase (afternoon), the coefficient is smaller (e.g., 0.8) to maintain stability.

[0632] The forgetting gate also performs rhythm-aware adjustment, adaptively adjusting the forgetting rate based on the IMF cycle length. The shorter the cycle, the faster the forgetting. Specific rules:

[0633] For cycles < 1 hour: basal forgetting rate × 0.3; for cycles 1 hour ≤ cycle < 6 hours: basal forgetting rate × 0.6; for cycles 6 hours ≤ cycle < 48 hours: basal forgetting rate × 0.85; for cycles ≥ 48 hours: basal forgetting rate × 0.95; rhythm stability is also considered. Highly stable rhythms (e.g., diurnal rhythms, stability > 0.9) reduce the forgetting rate by 10%; unstable rhythms increase the forgetting rate by 10%. The final forgetting threshold is: ft = basal value × cycle coefficient × stability coefficient.

[0634] Candidate memories are no longer a single value, but are decomposed into three components:

[0635] Trend component: Captures slowly changing baselines, calculated using IMF7-8 and residuals.

[0636]

[0637] Rhythmic components: capture periodic variations, calculated using IMF1-6.

[0638] For each IMFk: Where Ak is the amplitude weight.

[0639] Coupled components: capturing the interactions between indicators, based on coupled network computation.

[0640]

[0641] Final candidate memory:

[0642] This decomposition ensures that different types of information are processed independently, avoiding mutual interference.

[0643] An attention mechanism, or attention layer, is added on top of the LSTM layer to handle the delayed causal relationship between physiological indicators.

[0644] For each historical time step t', extract the hidden state vector ht' (from the third LSTM layer, 128-dimensional). Use the current hidden state ht to obtain the query vector through a linear transformation: Q = WQ × ht, where WQ is a 128 × 64 weight matrix and Q is 64-dimensional. For each historical time step: Kt' = WK × ht', where WK is also 128 × 64. Dot product attention: score base (t') = Q·Kt' / √64. Dividing by √64 is a scaling factor to prevent the gradient from vanishing due to an excessively large dot product. Causal delay modulation is added. For each possible source index j: obtain its Granger causality index GCIj and optimal delayj; calculate time alignment: align = exp(-|prediction time - (t' + delayj)| / σ); delay modulation score: score delay (t',j) = GCIj × align; Comprehensive attention score; score(t') = score base (t')+Σjscore delay (t',j); Normalization. Use the softmax function: attention(t') = exp(score(t')) / Σt”exp(score(t”)).

[0645] The attention layer consists of four attention heads: immediate attention head, short-term attention head, medium-term attention head, and long-term attention head. Each attention head has independent parameters and focus points.

[0646] Real-time attention head: query, key, value matrix dimension: 128×16 (generating a 16-dimensional representation), attention window: ±5 minutes, σ value: 0.5 (high concentration), primarily capturing real-time effects such as breathing.

[0647] Mid-term attention head: Matrix dimensions: 128×16, attention window: ±2 hours, σ value: 1.0, capturing short-term effects such as exercise and eating.

[0648] Mid-term attention head: Matrix dimensions: 128×16, attention window: ±8 hours, σ value: 2.0, capturing mid-term effects such as body temperature and hormones.

[0649] Long-term attention head: Matrix dimensions: 128×16, attention window: ±24 hours, σ value: 4.0, capturing long-term effects such as sleep accumulation.

[0650] The outputs of the four attention heads are concatenated into a 64-dimensional vector (4×16). This vector is then mapped back to 128 dimensions using the output matrix WO (64×128): Output = WO × Concat[head1,head2,head3,head4];

[0651] Add residual connections: Final output =Output + Primitive Hidden State

[0652] Layer normalization: for Final output Normalize it so that its mean is 0 and its variance is 1.

[0653] The decoder is designed to integrate all information and generate the final prediction. It integrates the outputs of four LSTM layers and an attention layer: LSTM layer 1 output: 32 dimensions; LSTM layer 2 output: 64 dimensions; LSTM layer 3 output: 128 dimensions; LSTM layer 4 output: 64 dimensions; attention layer output: 128 dimensions; current coupling effect: 150 dimensions (5×5 index pairs × 6 delay); total: 566 dimensions.

[0654] Use a gating mechanism to control the contributions of each part:

[0655] Gating value calculation: g = sigmoid(Wg × [all inputs] + bg), where Wg is a 566 × 5 matrix, generating 5 gating values.

[0656] Weighted fusion: Fusion = g1×LSTM1 + g2×LSTM2 + g3×LSTM3 + g4×LSTM4 + g5×Attention;

[0657] This gating mechanism allows the network to automatically learn the importance of different information sources.

[0658] The decoding layer consists of a first fully connected layer and a second fully connected layer.

[0659] First fully connected layer: Input: fused 566-dimensional vector, Output: 256-dimensional, Weight matrix: 566×256=145,096 parameters, Activation function: ReLU(x)=max(0,x), Dropout: randomly deactivates 20% of neurons during training.

[0660] The second fully connected layer: Input: 256-dimensional, Output: 128-dimensional, Weight matrix: 256×128=32,768 parameters, Activation function: ReLU, Dropout: 10%.

[0661] The output layer is divided into two branches:

[0662] Prediction branch: 8 output heads, each predicting the future value of an IMF. Each head is a 128×1 linear layer with no activation function (positive and negative values ​​are allowed). Post-processing: The 8 IMFs are summed and the residual is added to obtain the final predicted value.

[0663] Uncertainty branch: outputs the standard deviation of the prediction, a 128×1 linear layer, Softplus activation: softplus(x)=log(1+exp(x)), ensuring the output is positive, used to calculate the confidence interval: 95% interval = predicted value ± 1.96×standard deviation.

[0664] The first long short-term memory network model configured as described above can adapt to the multi-scale rhythm decomposition results and physiological rhythm coupling network, resulting in more accurate prediction results that are closer to user needs.

[0665] Initialization of the feature extraction layer (LSTM layer): The weight matrix is ​​initialized using Xavier, sampled from a uniform distribution U(-√(6 / (input dimension + output dimension)),√(6 / (input dimension + output dimension))). This ensures that the variance of the gradients remains consistent during forward and backward propagation. The forget gate bias is initialized to 1.0: This makes the network initially inclined to remember information rather than forget it, which helps in learning long-term dependencies. The biases of other gates are initialized to 0.

[0666] Attention layer initialization: The query, key, and value matrices are initialized using a scaled Xavier algorithm with a scaling factor of 1 / √number of heads = 0.5. This prevents the output of multi-head attention from becoming too large. The output matrix is ​​initialized using a normal distribution: N(0, 0.02). Smaller initial values ​​ensure the dominance of residual connections.

[0667] Fully connected layer initialization: He initialization is used: sampling from N(0, √(2 / input dimension)). This is suitable for the ReLU activation function and prevents gradient vanishing. The bias is initialized to 0.01 (a small positive value) to avoid dead neurons.

[0668] Create a complete initial long short-term memory network model, and assemble the complete model according to the following steps:

[0669] Create a model container (using PyTorch's nn.Module or TensorFlow's Model class). Add the layers in sequence: input preprocessing layer (normalization, encoding), four LSTM layers (parallel processing), attention layer, fusion layer, two fully connected layers, and output layer (two branches). Define the forward propagation function to specify how data flows through the network.

[0670] Total number of model parameters: LSTM layer: approximately 150,000 parameters; Attention layer: approximately 30,000 parameters; Fully connected layer: approximately 180,000 parameters; Total: approximately 360,000 parameters; Save the initial Long Short-Term Memory network model to a file as the starting point for training.

[0671] Data batch creation: Batch size: 32 samples; each sample: 48 hours input → 24 hours output; samples are generated using a sliding window, with the window moving 1 hour at a time; 10 days of data can generate approximately 200 training samples; random phase shift: ±2 hours, simulating biological clock shift; amplitude scaling: 0.8-1.2 times, simulating individual differences; noise addition: Gaussian noise, σ = 0.05 × signal amplitude; random masking: 10% of the time period is set to missing, for training robustness; batch randomization: the sample order is shuffled at the beginning of each epoch to prevent the model from memorizing the training order.

[0672] Phase 1: Layered pre-training. The training is divided into three phases, with the first phase being independent pre-training for each layer.

[0673] Pre-train the first LSTM layer (ultra-short rhythm hierarchical): Freeze all other layers and train only the first layer. Prepare dedicated training data containing only IMF1-2. Define a simplified loss function: L = MSE(predicted IMF1-2, actual IMF1-2), using the Adam optimizer with a learning rate of 0.001, β1 = 0.9, and β2 = 0.999. Train for 200 epochs, iterating through all data in each epoch. Monitor the validation set loss; if it doesn't decrease for 5 consecutive epochs, halve the learning rate. Save the trained parameters of the first layer.

[0674] Similarly, pre-training short-term rhythm layering, diurnal rhythm layering, and long-term rhythm layering is performed as follows: Second layer: 300 epochs, learning rate 0.001; Third layer: 500 epochs, learning rate 0.0005 (more important, more detailed training); Fourth layer: 200 epochs, learning rate 0.001.

[0675] Phase Two: Joint Fine-Tuning. After pre-training, perform end-to-end joint training. Load all pre-trained layer parameters. Unfreeze all layers to allow simultaneous updates. Use the full composite loss function:

[0676] L total =L mse +0.3×L rhythm +0.2×L phase +0.2×L couple +0.1×L uncertainty

[0677] Where: L mse : Mean squared error between predicted and actual values; L rhythmLoss of consistency in rhythmic characteristics (amplitude, period); L phase Phase continuity loss, penalizing phase jumps; L couple : Coupling consistency loss, ensuring adherence to causal relationships; L uncertainty Uncertainty calibration loss.

[0678] Use a small learning rate of 0.0001 to prevent corrupting pre-trained features. Apply gradient clipping, limiting the gradient norm to ≤5.0 to prevent gradient explosion. Train for 1000 epochs, periodically evaluating on the validation set. Use a learning rate schedule: 0.0001 for the first 300 epochs, 0.00005 for 300–700 epochs, and 0.00001 for 700–1000 epochs, preserving the best model (the one with the smallest loss on the validation set).

[0679] Phase 3: Personalized adaptation, using data from specific users for personalization.

[0680] Load the jointly trained model. Step 8.3.2: Freeze the bottom layer (the first three LSTM layers), and only update the top layer and decoder. Use the user's data from the last 7 days. Minimum learning rate of 0.00001, fine-tuning for 500 epochs. Focus on optimizing individual-specific parameters: baseline physiological values ​​(individual baselines such as heart rate and body temperature), rhythm amplitude (individual diurnal variation), phase preference (early bird or night owl pattern), and coupling strength (individual correlation strength between physiological indicators).

[0681] Starting with physiological state data from physiological state vectors, this study reveals multi-timescale patterns through multi-scale rhythm decomposition. It then uses a physiological rhythm coupling network to understand the interactions between physiological indicators. A specialized LSTM model is designed to integrate all information, and after training, an accurate physiological state prediction model is obtained. This physiological state prediction model is not a black box; each step has a clear physiological meaning, and the prediction results are interpretable and verifiable, providing precise prediction data for realizing biorhythm-based smart home environment regulation.

[0682] In one embodiment of this application, a PC algorithm based on the Pearl causal inference framework is used to discover the causal structure between environmental parameters and physiological state data, and to construct an environmental-physiological causal relationship model, based on physiological state data, environmental parameters, and a digital twin model of biological rhythms.

[0683] Initial relational variables are obtained based on physiological state data, environmental parameters, and the multi-scale biorhythm decomposition results in the biorhythm digital twin model.

[0684] Based on the initial relational variables and the biological rhythm numerical twin model, an initial complete graph is constructed, which includes multiple first-preset edges;

[0685] Independence tests were conducted on physiological state data based on the physiological rhythm coupling network and initial relational variables in the biorhythm digital twin model, and the first independence test results were obtained.

[0686] Based on the independence test results and the physiological rhythm coupling network in the biorhythm digital twin model, the multiple first preset edges of the initial complete graph are adjusted to obtain the skeleton graph, which includes multiple second preset edges.

[0687] Identify V-shaped triplets on the skeleton diagram to obtain V-shaped triplet data;

[0688] The orientation of some second preset edges in the skeleton diagram is identified by the direction propagation rule and the V-shaped structure triplet, and the orientation of the obtained second preset edges is determined as the first orientation.

[0689] Based on physiological state data and a digital twin model of biological rhythms, the direction of the second preset edge line of other parts is identified, and the direction of the second preset edge line of other parts is determined as the second direction;

[0690] Based on the skeleton diagram, the first direction, and the second direction, the causal structure is obtained;

[0691] Construct structural equations based on the causal structure;

[0692] The parameters of the structural equation model are evaluated based on physiological state data to obtain a parameter vector;

[0693] The standardized causal strength is obtained from the parameter vector;

[0694] Based on the multi-scale biorhythm decomposition results and standardized causal strength of the biorhythm digital twin model, the comprehensive scale causal strength is obtained;

[0695] Based on the causal structure and structural equations, an environmental-physiological causal relationship model is generated.

[0696] In this embodiment, to perform PC algorithm analysis, all data needs to be integrated into a unified data matrix. This data matrix D has a dimension of N×M, where N is the number of time points (e.g., 10080) and M is the number of variables. The variables include 11 directly observed variables (5 physiological indicators + 6 environmental parameters) and multi-scale rhythm decomposition results.

[0697] The specific data matrix construction process is as follows: Each physiological indicator is decomposed into rhythm components and residuals, using the IMF decomposition results of the multi-scale rhythm decomposition. For heart rate HR(t), based on the decomposition results of the multi-scale rhythm decomposition, it can be expressed as:

[0698] HR(t) = HR trend(t)+HR circadian (t)+HR ultradian (t)+HR residual (t)

[0699] Among them HR trend It is a long-term trend (corresponding to IMF7 and IMF8), HR circadian It is a circadian rhythm component (corresponding to IMF5 and IMF6), HR ultradian It is an ultrashort rhythm component (corresponding to IMF1 to IMF) 45 ), HR residual It is a random component.

[0700] In this way, the original heart rate time series is decomposed into four components, each serving as an independent variable in the PC algorithm. The same decomposition is performed on all five physiological indicators, resulting in 20 physiological rhythm variables. Adding six environmental variables, a total of 26 variables constitute the input to the PC algorithm.

[0701] The PC algorithm starts with a fully connected undirected graph, assuming that causal relationships may exist between all variables. In our scenario, this means constructing an initial complete graph with a 26×26 variable matrix and 325 possible edges (26×25 / 2). However, based on the physiological rhythm coupling network, we can intelligently reduce the initial complete graph. The Granger causality matrix (GC) from the Granger causality analysis results is shown below. matrix We have already identified which physiological indicators have significant coupling relationships. For example, if GC... matrix [HR][HRV] = 0.65, indicating that heart rate has a strong causal effect on heart rate variability; if GC matrix [BT][SBP] = 0.08, indicating that body temperature has a very weak effect on blood pressure.

[0702] Based on this physiological rhythm coupling network, we adopted a hierarchical strategy when constructing the initial graph:

[0703] First level: It is assumed that there is no direct causal relationship between environmental variables (environmental parameters are independently controlled);

[0704] Second layer: from environmental variables to physiological variables, retaining all possible edges (6×20=120);

[0705] The third layer: Among physiological variables, only the identified strong coupling relationships are retained (GC). matrix Edges with values ​​> 0.3.

[0706] In this way, the initial complete graph is reduced from 325 edges to about 150 edges, which greatly improves the efficiency of subsequent analysis.

[0707] The conditional independence test is the core of the PC algorithm, used to determine whether two variables are independent given other variables. For any two variables X and Y (which can be environmental or physiological variables), and a condition set S (a subset of other variables), the test checks whether X and Y are conditionally independent given S. The partial correlation coefficient is used as a measure of independence.

[0708] ρ(X,Y|S)=[ρ(X,Y)-ρ(X,S)×ρ(Y,S)] / [√(1-ρ 2 (X,S))×√(1-ρ 2 (Y,S))]

[0709] Where ρ(X,Y) is the Pearson correlation coefficient between X and Y, calculated using the following formula:

[0710]

[0711] Among them, X i X is the value of X at the i-th time point. - ρ(X,S) is the mean of X; ρ(X,S) is the multiple correlation coefficient between X and the condition set S; ρ(Y,S) is the multiple correlation coefficient between Y and the condition set S.

[0712] Taking temperature (TEMP) and heart rate (HR) as examples: Extract time series of TEMP(t) and HR(t) from the S1 data, with 10080 data points in each series. Calculate the direct correlation coefficient ρ(TEMP,HR) by multiplying the corresponding values ​​of the two time series, summing the results, and dividing by the product of their standard deviations. Assume a result of 0.45, indicating a moderate positive correlation. Based on the circadian rhythm coupling network, select possible confounding variables, such as body temperature (BT) (because temperature affects body temperature, and body temperature affects heart rate). Calculate the conditional correlation coefficient ρ(TEMP,HR|BT). First, calculate ρ(TEMP,BT) = 0.72 and ρ(HR,BT) = 0.58, then substitute them into the partial correlation formula:

[0713] ρ(TEMP,HR|BT)=(0.45-0.72×0.58) / [√(1-0.72 2 )×√(1-0.58 2 )]

[0714] =0.032 / (0.69×0.81) =0.057

[0715] Since the conditional correlation coefficient of 0.057 is less than the threshold of 0.1, temperature and heart rate are independent under a given body temperature, meaning that the effect of temperature on heart rate is mediated by body temperature.

[0716] The independence test is enhanced by using the results of multi-scale biological rhythm decomposition, which allows the independence test to be performed separately at different time scales, thus greatly improving the accuracy of causal discovery.

[0717] For the ultrashort rhythm component (5–20 minute cycles), a shorter time window (1 hour) is used for independence testing. This is because ultrashort rhythms reflect rapid physiological regulation, and a longer window would mask these rapid changes. Specifically, the 10080 data points are divided into 168 1-hour segments, the correlation coefficient is calculated for each segment, and then the average is taken. For the diurnal rhythm component (24-hour cycle), complete 7 days of data are used for independence testing. Simultaneously, considering the phase difference between diurnal and circadian rhythms, a time delay is introduced when calculating the correlation coefficient. If the user's body temperature rhythm phase is 2 hours ahead of the heart rate rhythm, the body temperature data is shifted forward by 2 hours when calculating ρ(BT_circadian, HR_circadian).

[0718] Multiscale independence tests can uncover different types of causal relationships. For example, the immediate effects of temperature on heart rate (through skin vascular responses) and its long-term effects (through thermoregulation) can be identified separately.

[0719] Based on the results of the conditional independence test, edges in the initial complete graph that do not represent a causal relationship are removed. Zero-order independence test (unconditional): This directly tests the correlation between each pair of variables. For example, testing the direct correlation coefficient between illumination (LUX) and blood pressure (SBP). If ρ(LUX,SBP) = 0.03 < 0.1, the edge LUX-SBP is directly removed. At this stage, approximately 30% of the edges can usually be removed, primarily those between environmental variables and irrelevant physiological indicators.

[0720] First-order independence test (conditional on a single variable): For each remaining edge, test whether there exists a third variable such that the two edge variables are conditionally independent. For example, for the TEMP-HR edge, test the conditional independence when conditional on each variable such as BT, HRV, SBP, and DBP. If ρ(TEMP,HR|BT) < 0.1, record BT as the separating set and delete the TEMP-HR edge.

[0721] Second-order independence test (conditional on two variables): For the edges that are still retained, test the independence conditional on the two variables. For example, test ρ(LUX,HR|CCT,BT), that is, whether illuminance and heart rate are independent given color temperature and body temperature. This stage is mainly used to identify complex confounding relationships.

[0722] Higher-order independence tests: Continue increasing the size of the condition set until no more edges can be deleted. In practice, testing up to order 3 or 4 is usually sufficient, as higher-order condition sets lead to a decrease in statistical power.

[0723] The Granger causality matrix from the Granger causality analysis results shows that variable A strongly influences variable B (GC). matrix If [A][B]>0.5), then when testing any edge related to B, we prioritize adding A to the condition set. For example, if body temperature strongly affects heart rate variability, when testing the TEMP-HRV edge, we first test ρ(TEMP,HRV|BT).

[0724] Delay matrix matrix It also provides important information. If Delay matrix [A][B] = 15 minutes, meaning that the effect of A on B has a 15-minute delay. When conducting conditional independence tests, we shift the data for A forward by 15 minutes before calculating the correlation coefficient, which allows us to more accurately capture causal relationships.

[0725] A V-structure (also known as a convergent structure) is a structure of the form A→C←B, where A and B are not directly connected. Identifying V-structures is crucial for determining causal directions. In the skeleton graph after edge removal, find all triples (A,C,B) where AC and BC are connected, but A and B are not. For each such triple, check whether A and B become dependent (rather than independent) given C.

[0726] Suppose we have a triple (TEMP, BT, HR), where TEMP is connected to body temperature (BT), and HR is also connected to body temperature (BT), but TEMP and HR are not directly connected. Calculate:

[0727] ρ(TEMP,HR) = 0.05 (almost independent);

[0728] ρ(TEMP,HR|BT) = 0.38 (becomes relevant given BT);

[0729] This phenomenon of "becoming dependent after being given an intermediate variable" strongly suggests the existence of a V-structure TEMP→BT←HR. This means that both temperature and heart rate affect body temperature, rather than body temperature affecting both temperature and heart rate simultaneously (which is also physically illogical).

[0730] If a physiological indicator has a particularly strong diurnal rhythm (such as body temperature, which accounts for 40% of the total variation), then this indicator is more likely to be a convergence point of the V-structure, because multiple factors jointly regulate strong rhythm indicators.

[0731] After identifying the V-structure, the first direction of the other edges is determined using direction propagation rules. The main rules include:

[0732] Rule 1: If there exists a directed path A→B→C and AC are connected, then the direction must be A→C (to avoid forming a new V-structure).

[0733] Rule 2: If there exists A→BC and A and C are not connected, then the direction must be B→C (to avoid forming a new V-structure).

[0734] Rule 3: Directed cycles cannot be formed (causality cannot be cyclical).

[0735] When applying these rules, physiological rhythm coupled networks should be given priority. If the Granger causality matrix GC... matrix [A][B] is significantly greater than GC matrix [B][A], prioritize orienting the edge to A→B.

[0736] For edges whose direction cannot be determined by the PC algorithm, temporal information and domain knowledge are used for orientation to obtain a second direction.

[0737] Time priority principle: Events that occur earlier in time are more likely to be the cause. Calculate the cross-correlation function:

[0738] R XY (τ)=Σ t [X(t)×Y(t+τ)] / [σ X ×σ Y ×N]

[0739] Where τ is the time delay, σ is the standard deviation, and N is the number of data points. If R XY (τ) reaches its maximum value when τ>0, indicating that X leads Y, and it is more likely that X→Y.

[0740] For example, calculating the cross-correlation function between temperature and body temperature reveals Rw TEMP The BT value is the highest at 25 minutes, indicating that temperature changes have the strongest impact on body temperature after 25 minutes. Therefore, the direction is determined to be TEMP→BT.

[0741] Delay matrix matrix It directly provides causal delay information. If Delay matrix [HR][HRV] = 5 minutes, while Delay matrix If [HRV][HR] = 0 (indicating no reverse effect), then the direction is determined to be HR→HRV.

[0742] If the circadian rhythm phase of variable A leads that of variable B, and their circadian rhythms are correlated, then A is more likely to influence B. For example, if the circadian rhythm phase of body temperature is 2 hours earlier than that of heart rate, and their diurnal components are highly correlated (ρ>0.7), then the direction of BT→HR is supported.

[0743] Based on the skeleton diagram, the first direction, and the second direction, the causal structure is obtained;

[0744] After determining the causal structure, it is necessary to quantify the strength of each causal edge. Using structural equation modeling, the causal effect of each edge A→B is estimated.

[0745] Model form: B(t)=β0+β1×A(t-δ)+Σ i γ i ×C i (t)+ε(t)

[0746] Where: A(t-δ) is the causal variable, δ is the causal delay; C i These are the other parent nodes of B (nodes pointing to B in the causal graph); β1 is the direct causal effect coefficient of A on B; γ i ε is the effect coefficient of other parent nodes; ε is the error term.

[0747] Using physiological state data, the parameter is estimated using the least squares method: β=(X'X) -1 X'Y

[0748] Where X is the design matrix (containing A and all C) i Y is the data vector of B.

[0749] Causal strength is defined as the standardized effect coefficient: Causal Strength =|β1|×(σ A / σ B )

[0750] This standardization ensures that the causal strength is between 0 and 1, making it easier to compare the causal strength of different pairs of variables.

[0751] For each causal edge A→B, calculate: Short-term causal strength: calculate the causal effect using the IMF1-IMF4 components of A and B; Diurnal causal strength: calculate the causal effect using the IMF5-IMF6 components of A and B; Long-term trend causal strength: calculate the causal effect using the IMF1-IMF4 components of A and B. i -IMF8 component calculation of causal effects.

[0752] Integrating all the analysis results to form the final model, after the PC algorithm is completed, a complete environmental-physiological causal relationship model is obtained, which includes: a causal graph structure G = (V, E), where V is 26 nodes (6 environmental + 20 physiological rhythm variables) and E is a set of definite causal edges.

[0753] For each edge e = (i,j) ∈ E, the causal parameter matrix P stores the causal strength S. ij (between 0 and 1), causal delay D ij (minutes), confidence level Cij (Based on statistical significance), dynamic causal equations, which describe the dynamics of the entire system:

[0754] dY i / dt=Σ j S ji ×X_j(tD ji )+Σ k G ki ×Y k (t)+R i (t)

[0755] Where Y i X is the i-th physiological variable. j It's an environment variable, G ki It is the physiological rhythm coupling network provided by S2, R i It is a rhythm driver.

[0756] In one embodiment of this application, the multi-objective optimization problem includes: a health promotion objective function, a comfort objective function, an energy consumption cost objective function, and an equipment lifespan objective function; the multi-objective optimization problem is transformed into a quadratic unconstrained binary optimization function, and the Pareto optimal solution set is obtained by applying a quantum optimization algorithm to the quadratic unconstrained binary optimization function, including:

[0757] Based on the multi-objective optimization problem, the binary optimization matrices are obtained;

[0758] The constraints are transformed into penalty terms and added to each binary optimization matrix to obtain a quadratic unconstrained binary optimization function;

[0759] By mapping the quadratic unconstrained binary optimization function to the Ising model and gradually reducing the transverse field intensity of the Ising model until it finally stabilizes at the lowest energy state, the Pareto optimal solution set is obtained.

[0760] In the embodiments of this application, the multi-objective optimization problem (QUBO problem) uses binary variables (0 and 1), while quantum physical systems naturally use spin variables (upward +1 and downward -1). The Ising model is the standard physical model describing spin systems, and quantum annealing hardware is a physical system that simulates the Ising model. Therefore, it is necessary to establish a mathematical equivalence between QUBO and the Ising model.

[0761] The standard form of the QUBO problem is: minimize f(q) = Σ ij Q ij q i q j ;

[0762] Where q is a binary vector, each q i It takes the value 0 or 1, and Q is the coefficient matrix.

[0763] The energy function of the Ising model is: E(s) = -Σ j J j s i s j -Σ i h i s i ;

[0764] Where s is the spin vector, and each s i The value can be +1 or -1, where J is the coupling strength matrix and h is the local field vector.

[0765] Mapping is achieved through variable substitution: q i =(s i +1) / 2

[0766] The meaning of this substitution is: when the spin is upward s i When =+1, it corresponds to binary q i =1; when spin down s i When = -1, it corresponds to binary q. i =0; Substitute the substitution into the QUBO function, and after algebraic expansion:

[0767] f(q)=Σ j Q j ×((s i +1) / 2)×((s j +1) / 2)=(1 / 4)Σ ij Q ij (s i s j +s i +s j +1)

[0768] Reorganizing the terms, we obtain the Ising form: f(s) = -Σ j J j s i s j -Σ i h i s i +constant

[0769] The mapping relationship is: J ij =-Q ij / 4 (coupling term); h i =-(Σ j Q ij ) / 2 (local field);

[0770] The transverse field is the core concept of quantum annealing. In classical physics, spin can only go up or down. However, in quantum physics, spin can exist in a superposition state, possessing both up and down probabilities. The transverse field is the "driving force" that causes spin to undergo quantum superposition. The stronger the transverse field, the greater the degree of superposition, and the easier it is to transition between different states. The complete process of quantum annealing is as follows:

[0771] Initial stage (t=0): The transverse field strength Γ0 is set to its maximum value (e.g., 10), all spins are in a complete superposition state, the system energy is dominated by the transverse field, and the problem energy (Ising energy) has little influence. Physical picture: the system is in a "quantum liquid" state, and all possible solutions are explored equally. Intermediate stage (t=0) <t<T anneal As the transverse field strength gradually decreases: Γ(t)=Γ0×(1-t / T_anneal), the energy of the problem gradually becomes important. It begins to "sensor" which states have lower energies, and the physical picture is that the system gradually "crystallizes" from a "liquid" state, and begins to tend towards a lower energy configuration.

[0772] During the annealing process, quantum tunneling enables the system to cross the energy barrier, which is the key to the superiority of quantum optimization over classical optimization.

[0773] Quantum annealing: It can directly "pass through" the energy barrier (quantum tunneling), making it easier to find the global optimum. Specific mechanism:

[0774] Suppose the system is currently in state A (local optimum) and the global optimum is in state B, with an energy barrier in between.

[0775] Quantum method: tunneling directly from A to B without passing through a high-energy intermediate state.

[0776] Final stage (t=T) anneal ); The transverse field decreases to zero: Γ(T) anneal When ) = 0, quantum superposition disappears, and the system "collapses" to a definite classical state, which corresponds to the ground state (lowest energy state) of the Ising model. The optimal solution of QUBO is obtained through inverse mapping.

[0777] Measurement and decoding: Measure the direction of each spin (+1 or -1), via q i =(s i +1) / 2 converts back to binary and decodes the binary sequence into environmental control parameters.

[0778] A single quantum annealing operation yields only one solution, while multi-objective optimization problems require a set of Pareto optimal solutions to choose from. By changing the combination of objective functions and running quantum annealing multiple times, a solution set can be generated.

[0779] Select representative points in the four-dimensional weight space:

[0780] w1 (health weight): from 0.1 to 0.7, step size 0.1; w2 (comfort weight): from 0.1 to 0.7, step size 0.1; w3 (energy consumption weight): from 0.1 to 0.7, step size 0.1; w4 (device weight): 1-w1-w2-w3 (guaranteed sum to 1); This can generate approximately 50 to 100 different weight combinations. For each weight combination (w1, w2, w3, w4): construct a combination QUBO matrix:

[0781] Q combined = w1×Q1+w2×Q2+w3×Q3+w4×Q4

[0782] Q1 is the QUBO matrix for health goals, Q2 is for comfort goals, and so on.

[0783] Mapping to the Ising model, using the aforementioned mapping formula, Q... combined Convert to Ising parameters J and h.

[0784] Perform quantum annealing, initialize the transverse field Γ0 = 10, run the annealing process for 1 second, gradually reduce the transverse field to 0, obtain the spin configuration s*, convert s to binary q, decode q into the environmental parameter sequence X, and calculate the performance of this scheme on four targets [f1*, f2*, f3*, f4*].

[0785] The quantum annealing algorithm described above successfully transforms complex multi-objective optimization problems into forms that can be efficiently solved by quantum computing, and generates diverse Pareto optimal solution sets, providing users with scientific and flexible environmental control solutions.

[0786] Secondly, embodiments of this application provide a control system for an environment and equipment cluster, including:

[0787] The first acquisition module is used to acquire physiological state data and environmental parameters, and to obtain a physiological state vector based on the physiological state data.

[0788] The first building module is used to perform empirical mode decomposition on physiological state vectors and construct a digital twin model of biological rhythms;

[0789] The second building module is used to discover the causal structure between environmental parameters and physiological state data based on physiological state data, environmental parameters and biological rhythm digital twin models, and to build an environmental-physiological causal relationship model.

[0790] The third construction module is used to construct a multi-objective optimization problem based on the biological rhythm digital twin model and the environmental physiological causal relationship model;

[0791] The first module is used to transform the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and to obtain the Pareto optimal solution set by applying the quantum optimization algorithm to the quadratic unconstrained binary optimization function.

[0792] The first regulation module is used to regulate the environment and equipment clusters based on the Pareto optimal solution set and the digital twin model of biological rhythms.

[0793] The functions of each module in each device in the embodiments of this application can be found in the corresponding descriptions in the above methods, and will not be repeated here.

[0794] Figure 2 A structural block diagram of an electronic device according to an embodiment of this application is shown. Figure 2 As shown, the electronic device includes a memory 410 and a processor 420. The memory 410 stores instructions that can be executed on the processor 420. When the processor 420 executes the instructions, it implements the environmental and device cluster control method described in the above embodiments. The number of memories 410 and processors 420 can be one or more. This electronic device is intended to represent various forms of digital computers, such as laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device can also represent various forms of mobile devices, such as personal digital processors, cellular phones, smartphones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the present application described and / or claimed herein.

[0795] The electronic device may also include a communication interface 430 for communicating with external devices and exchanging data. The devices are interconnected using different buses and can be mounted on a common motherboard or otherwise as needed. The processor 420 can process instructions executed within the electronic device, including instructions stored in or on memory to display graphical information of a GUI on an external input / output device (such as a display device coupled to the interface). In other embodiments, multiple processors and / or multiple buses can be used with multiple memories and multiple memory modules, if desired. Similarly, multiple electronic devices can be connected, each providing some of the necessary operations (e.g., as a server array, a group of blade servers, or a multiprocessor system). The bus can be divided into address buses, data buses, control buses, etc. For ease of illustration, Figure 2 The bus is represented by a single thick line, but this does not mean that there is only one bus or one type of bus.

[0796] Optionally, in a specific implementation, if the memory 410, processor 420 and communication interface 430 are integrated on a single chip, the memory 410, processor 420 and communication interface 430 can communicate with each other through an internal interface.

[0797] It should be understood that the aforementioned processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. General-purpose processors can be microprocessors or any conventional processor. It is worth noting that the processor can be a processor supporting the Advanced Reduced Instruction Set Computing (RISC) machine architecture (ARM).

[0798] This application provides a computer-readable storage medium (such as the memory 410 described above) that stores computer instructions, which, when executed by a processor, implement the method provided in this application.

[0799] Optionally, memory 410 may include a program storage area and a data storage area, wherein the program storage area may store the operating system and applications required for at least one function; the data storage area may store data created based on the use of the electronic device, etc. Furthermore, memory 410 may include high-speed random access memory, and may also include non-transitory memory, such as at least one disk storage device, flash memory device, or other non-transitory solid-state storage device. In some embodiments, memory 410 may optionally include memory remotely located relative to processor 420, and these remote memories can be connected to the electronic device via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0800] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this application. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of those different embodiments or examples.

[0801] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the order or number of the indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this application, "a plurality of" means two or more, unless otherwise explicitly specified.

[0802] The above are merely specific embodiments 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 variations or substitutions within the technical scope disclosed in this application, and these should all be included 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 controlling an environment and equipment cluster, characterized in that, The method comprises the following steps: acquiring physiological state data and environmental parameters, and obtaining a physiological state vector based on the physiological state data; performing empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model; applying a PC algorithm of a Pearl causal inference framework to discover a causal structure between the environmental parameters and the physiological state data based on the physiological state data, the environmental parameters, and the biological rhythm digital twin model, and constructing an environmental-physiological causal relationship model; constructing a multi-objective optimization problem according to the biological rhythm digital twin model and the environmental-physiological causal relationship model; transforming the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and obtaining a Pareto optimal solution set by performing quantum optimization algorithm on the quadratic unconstrained binary optimization function; based on the Pareto optimal solution set and the biological rhythm digital twin model, regulating the environment and the device cluster.

2. The method of claim 1, wherein, The biological rhythm digital twin model comprises a multi-scale biological rhythm decomposition result, a physiological rhythm coupling network, and a physiological state prediction model, and the method of performing empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model comprises the following steps: performing empirical mode decomposition on the physiological state vector to obtain a multi-scale biological rhythm decomposition result; performing Granger causality analysis on the multi-scale biological rhythm decomposition result to obtain a Granger causality analysis result; performing transfer entropy calculation on the multi-scale biological rhythm decomposition result to obtain a transfer entropy result; determining an influence type result according to the multi-scale biological rhythm decomposition result; constructing a physiological rhythm coupling network according to the Granger causality analysis result, the transfer entropy result, and the influence type result; obtaining a physiological state prediction model according to the multi-scale rhythm decomposition result, the physiological rhythm coupling network, and the physiological state vector.

3. The method of claim 2, wherein, The method of performing empirical mode decomposition on the physiological state vector to obtain a multi-scale biological rhythm decomposition result comprises the following steps: identifying the physiological state vector to obtain all local maximum points and all local minimum points; connecting the all local maximum points into a smooth upper envelope line and connecting the all local minimum points into a smooth lower envelope line by a cubic spline interpolation method, the upper envelope line and the lower envelope line surrounding the physiological state vector; calculating the average values of the upper envelope line and the lower envelope line to obtain a local mean curve; subtracting the local mean curve from the physiological state vector to obtain a first candidate intrinsic mode component; judging whether the first candidate intrinsic mode component satisfies an intrinsic mode condition; in the case that the first candidate intrinsic mode component satisfies the intrinsic mode condition, determining a first intrinsic mode; in the case that the first candidate intrinsic mode component does not satisfy the intrinsic mode condition, taking the first candidate intrinsic mode component as a physiological state vector, repeating the above process until the first candidate intrinsic mode component obtained satisfies the intrinsic mode condition, and obtaining a first intrinsic mode; subtracting the first intrinsic mode from the physiological state vector to obtain a first residual; The first residual is taken as a physiological state vector, and the above process is repeated until the obtained first residual becomes a monotonic function or the extreme points thereof are less than 2, and the first eigenmode in the repeated process is sequentially determined as a second eigenmode, a third eigenmode, a fourth eigenmode, a fifth eigenmode, a sixth eigenmode, a seventh eigenmode, and an eighth eigenmode; According to the first eigenmode, the second eigenmode, the third eigenmode, the fourth eigenmode, the fifth eigenmode, the sixth eigenmode, the seventh eigenmode, and the eighth eigenmode, a multi-scale biological rhythm decomposition result is obtained.

4. The method of claim 3, wherein, The analysis of the multi-scale biological rhythm decomposition result by Granger causality analysis includes: The multi-scale biological rhythm decomposition result is paired by Granger causality analysis based on the principle of similar periods, and the causal relationship between each biological rhythm is obtained; According to the multi-scale biological rhythm decomposition result, a first autoregressive model is constructed; The causal relationship between the corresponding biological rhythms is added to the first autoregressive model to determine a second autoregressive model; According to the first autoregressive model and the second autoregressive model, the causal strength is determined; According to the causal strength, the delay time is determined; According to the causal strength and the delay time, the Granger causality analysis result is obtained.

5. The method of claim 4, wherein, The physiological state prediction model is obtained according to the multi-scale rhythm decomposition result, the physiological rhythm coupling network, and the physiological state vector, which includes: According to the multi-scale rhythm decomposition result and the physiological rhythm coupling network, a first long short-term memory network model is constructed, the network model includes an input layer, a feature extraction layer, an attention layer, a fusion layer, and a decoding layer, the input layer includes an original value channel, a rhythm decomposition channel, a feature parameter channel, a coupling influence channel, and a time encoding channel, the feature extraction layer includes an ultra-short rhythm layer, a short-term rhythm layer, a circadian rhythm layer, and a long-term rhythm layer, the attention layer includes an instant attention head, a short-term attention head, a medium-term attention head, and a long-term attention head, and the decoding layer includes a first full connection layer and a second full connection layer; The first long short-term memory network model is initialized to obtain an initial long short-term memory network model; According to the physiological state vector, the initial long short-term memory network model is trained and tested to obtain a physiological state prediction model.

6. The method of claim 5, wherein, The PC algorithm of the Pearl causal inference framework is used to discover the causal structure between the environmental parameters and the physiological state data based on the physiological state data, the environmental parameters, and the biological rhythm digital twin model, and an environmental physiological causal relationship model is constructed, which includes: According to the physiological state data, the environmental parameters, and the multi-scale biological rhythm decomposition result in the biological rhythm digital twin model, variable data is obtained and a variable matrix is constructed; According to the variable matrix, an initial complete graph is constructed, and adjacent two variables in the initial complete graph have a first preset edge line connected between them; Performing an independence test on the physiological state data based on the physiological rhythm coupling network in the biological rhythm digital twin model and the initial relationship variable to obtain a first independence test result; Adjusting a plurality of first preset edges of the initial complete graph based on the independence test result and the physiological rhythm coupling network in the biological rhythm digital twin model to obtain a skeleton graph, wherein the skeleton graph comprises a plurality of second preset edges; Identifying a V-shaped structure triple in the skeleton graph to obtain V-shaped structure triple data; Identifying the direction of part of the second preset edges in the skeleton graph through a direction propagation rule and the V-shaped structure triple, and determining the direction of the part of the second preset edges as a first direction; Identifying the direction of other parts of the second preset edges based on the physiological state data and the biological rhythm digital twin model, and determining the direction of the other parts of the second preset edges as a second direction; Obtaining a causal structure based on the skeleton graph, the first direction, and the second direction; Constructing a structural equation based on the causal structure; Evaluating the parameters of the structural equation based on the physiological state data to obtain a parameter vector; Obtaining a standardized causal strength based on the parameter vector; Obtaining a comprehensive scale causal strength based on the multi-scale biological rhythm decomposition result of the biological rhythm digital twin model and the standardized causal strength; Generating an environmental physiological causal relationship model based on the causal structure and the structural equation.

7. The method of claim 6, wherein, The multi-objective optimization problem includes a health promotion objective function, a comfort objective function, an energy consumption cost objective function, and a device life objective function; the multi-objective optimization problem is converted into a quadratic unconstrained binary optimization function, and a set of Pareto optimal solutions is obtained by optimizing the quadratic unconstrained binary optimization function according to a quantum optimization algorithm, and the set of Pareto optimal solutions includes: Obtaining each binary optimization matrix according to the multi-objective optimization problem; Adding a constraint condition converted into a penalty term to each binary optimization matrix to obtain a quadratic unconstrained binary optimization function; Mapping the quadratic unconstrained binary optimization function to an Ising model, and gradually reducing the transverse field strength of the Ising model until it is finally stabilized in a state with the lowest energy to obtain a set of Pareto optimal solutions.

8. A control system for an environment and equipment cluster, characterized in that, The method comprises: A first acquisition module is configured to acquire physiological state data and environmental parameters, and obtain a physiological state vector based on the physiological state data; A first construction module is configured to perform empirical mode decomposition on the physiological state vector to construct a biological rhythm digital twin model; A second construction module is configured to discover a causal structure between the environmental parameters and the physiological state data by using a PC algorithm of a Pearl causal inference framework based on the physiological state data, the environmental parameters, and the biological rhythm digital twin model, and construct an environmental physiological causal relationship model; A third construction module is configured to construct a multi-objective optimization problem based on the biological rhythm digital twin model and the environmental physiological causal relationship model; The first obtaining module is configured to convert the multi-objective optimization problem into a quadratic unconstrained binary optimization function, and obtain a Pareto optimal solution set by performing quantum optimization algorithm on the quadratic unconstrained binary optimization function. The first regulation module is configured to regulate an environment and a device cluster based on the Pareto optimal solution set and the biological rhythm digital twin model.

9. An electronic device, comprising: The apparatus comprises: at least one processor; and a memory connected to the at least one processor in communication; wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the method of any one of claims 1-7.

10. A computer readable storage medium, the computer readable storage medium storing computer instructions, the computer instructions being executed by a processor to implement the method of any one of claims 1-7.

Images