An energy balance-based simulation modeling method for an asymmetric temperature control system

Abstract

The application discloses an energy balance-based asymmetric temperature control system simulation modeling method, and belongs to the technical field of industrial temperature control system modeling and simulation. Firstly, a temperature rising step experiment and a temperature falling step experiment are performed, steady-state data is acquired, a balance point power function is derived, equivalent verification is performed, a temperature rising curve and a temperature falling curve are derived, and identification parameters are derived. A net power deviation is defined, the identification parameters are subjected to smooth transition, gain parameters are acquired, a first-order linear differential model is established, a first-order backward difference method is adopted for discretization, recursion is performed based on a set simulation step, and a recursion model is acquired. A plurality of recursion models are derived based on repetition of the simulation step, a temperature response curve of the temperature control system is established, and simulation modeling is completed. The method only needs two step experiments, real-time parameters are calculated through the net power deviation, real smooth and continuous dynamic simulation is realized, and a reliable simulation method is provided for temperature digital twinning and temperature controller design.

Description

Technical Field

[0001] This invention relates to the field of industrial temperature control system modeling and simulation technology, and in particular to a simulation modeling method for asymmetric temperature control systems based on energy balance. Background Technology

[0002] In industrial temperature control systems, the asymmetry of dynamic characteristics between heating and cooling (or cooling) processes is a common phenomenon. Whether it's a heat-insulating device relying solely on natural heat dissipation with a heater (fast heating, slow cooling), a refrigeration device equipped with a compressor (fast cooling, slow heating), or even a system with bidirectional active regulation capabilities, the dynamics of heating and cooling often differ significantly. Therefore, accurate modeling and dynamic simulation of temperature control systems are crucial for controller design, parameter tuning, and the construction of digital twin platforms. Currently, existing modeling and simulation methods mainly suffer from the following shortcomings: 1. Pure energy balance mechanism models rely on physical parameters and have poor engineering applicability; For a natural heat dissipation temperature control system equipped only with a heater, its energy balance can be described as a dynamic balance between thermal inertia, heating power, and heat dissipation power. While the model can accurately reflect the system's thermal dynamics, it requires prior acquisition or identification of the total heat capacity C and equivalent thermal resistance R. th Heating efficiency η h These are internal physical parameters. These parameters are difficult to measure directly in actual engineering projects, making model identification complex, limiting their applicability, and hindering their widespread adoption.

[0003] 2. A single FOPDT model cannot simultaneously adapt to the asymmetric dynamics of heating and cooling; Existing technologies often use a first-order inertial plus pure time delay (FOPDT) model to describe the dynamics of temperature-controlled objects, identifying a set of model parameters through a single step test. However, the time constants τ and pure time delay L in the heating and cooling directions often differ significantly. A single FOPDT model can only match the dynamics in a single direction and cannot accurately reproduce the characteristics of the entire heating and cooling process, making it difficult to meet the requirements of high-precision simulation.

[0004] 3. Hard switching between the two models is not smooth and lacks a switching benchmark; Existing technologies employ a segmented approach, constructing and switching heating and cooling models independently, in an attempt to adapt to asymmetric characteristics. However, this method lacks a unified model switching benchmark, with switching thresholds often relying on empirical settings and failing to use net power deviation as the basis for parameter switching. Near the switching point, issues such as abrupt temperature changes and discontinuous simulation curves easily arise. It cannot realistically reflect the smooth transition characteristics of the physical system, making it difficult to support reliable digital twin simulation and controller development because the model does not meet the requirements of real-world operating conditions.

[0005] Therefore, existing technologies still need to be improved and enhanced. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the purpose of this invention is to provide a simulation modeling method for asymmetric temperature control systems based on energy balance. This method aims to solve the problems of existing technologies that often use first-order inertial plus pure time delay (FOPDT) models to describe the dynamics of temperature-controlled objects, which can only match dynamics in a single direction and cannot accurately reproduce the characteristics of the entire dynamic process of heating and cooling. Furthermore, the dual-model approach lacks a unified model switching benchmark, and the switching threshold often depends on empirical settings and does not use net power deviation as the basis for parameter switching, resulting in simulation models that do not conform to actual working conditions.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A simulation modeling method for an asymmetric temperature control system based on energy balance, the simulation modeling method includes: Step 1: Based on the physical system setting of heating power, conduct heating step experiment and cooling step experiment to obtain steady-state temperature and power duty cycle data, and obtain the equilibrium point power function; Step 2: Perform equivalent verification of the equilibrium point power function based on the first-order inertial plus pure time delay model, and obtain the heating curve and cooling curve based on the equilibrium point power function to obtain the identification parameters. Step 3: Define the net power deviation, and based on the defined net power deviation, smooth the transition of the identification parameters to obtain the gain parameters. Establish a first-order linear differential model according to the equilibrium point power function, discretize it using the first-order backward difference method, and recursively deduce based on the set simulation step size to obtain the recursive model. Step 4: Repeat step 3 based on the simulation step size to obtain multiple sets of recursive models, establish the temperature response curve of the temperature control system, and complete the simulation modeling.

[0008] Furthermore, in step 1, the temperature step experiment includes: Set a small heating power in the physical system, for example, u=10%, and record the initial temperature T after steady-state heating. h0 Initial power duty cycle u h0 Initial ambient temperature T amb_h0 Then, the heating power is switched to a higher level, for example, u=100%, and the temperature response curve T of the entire heating process is recorded. h (t), until the system stabilizes again, and record the final steady-state temperature T. h1 Final heating power duty cycle u h1 , steady-state ambient temperature T amb_h1 ; The cooling step test includes: Set a relatively large heating power in the physical system, for example, u=100%, and record the initial cooling temperature T after steady-state cooling. c0Initial power duty cycle for cooling, u c0 Initial ambient temperature T during cooling amb_c0 Then, the heating power is switched to a lower level, for example, u=10%, and the temperature response curve T of the entire cooling process is recorded. c (t), until the system stabilizes again, and record the final cooling steady-state temperature T. c1 Final cooling power duty cycle u c1 Cooling steady-state ambient temperature T amb_c1 .

[0009] Furthermore, in step 1, based on the steady-state temperature and power duty cycle data, the equilibrium point power function is obtained, and the energy balance equation is first established: C·(dT / dt)=η h ·P max ·u-(TT amb ) / R th ; Where T is the system temperature, t is time, C is the total heat capacity of the system, and η is the system temperature. h For heating efficiency, P max T represents the heater's maximum power, u is the power duty cycle (the ratio of actual heating power to maximum power), and its value ranges from 0 to 1 or from 0% to 100%. amb R represents the ambient temperature. th The equivalent thermal resistance between the system and the environment; In steady state, dT / dt = 0. Substituting this into the energy balance equation, we obtain the equilibrium power u. bal The formula is as follows: u bal =1 / (η h ·P max ·R th )·TT amb / (η h ·P max ·R th ); Let α = 1 / (η) h ·P max ·R th If the equilibrium point power u is... bal Represented as: u bal =αT-αT amb ; When the ambient temperature is constant, the equilibrium point power function is expressed as a univariate linear equation as follows: u bal (T) = αT + γ; When the ambient temperature changes, the equilibrium power function is expressed as a two-variable linear equation: u bal(T,T amb )=αT+βT amb .

[0010] Furthermore, the coefficients in the equilibrium point power function are obtained using either a directional calculation method or the least squares method. The directional calculation includes: Using the two steady-state points from the temperature increase step experiment in step 1, substitute them into the corresponding formal equation to obtain the coefficient (α) for the heating direction. h γ h ) or (α) h β h Similarly, the coefficient (α) of the cooling direction can be obtained by solving the two steady-state points of the cooling step test in step 1. c γ c ) or (α) c β c ); For the unary form, take α = (α h +α c ) / 2, γ=(γ h +γ c ) / 2; For the binary form, take α=(α h +α c ) / 2, β=(β h +β c ) / 2; Least squares methods include: Substitute all steady-state points into the corresponding function model and solve for the optimal coefficients directly using the least squares method.

[0011] Furthermore, in step 2, the equivalent verification of the equilibrium power function based on the first-order inertial plus pure time delay model includes: Considering a step input, at t=0, the heating power jumps from the initial power u0 to the final power u1 and remains thereafter; let the initial steady-state temperature be T0, satisfying u0=u bal (T0); Based on the equilibrium point power model from step 1, and using the initial conditions, the temperature response is derived as follows: T(t) = T0 + K·(u1-u0)·(1-e -(t-L) / τ ); Where τ=R th C is the time constant, K=η h P max R th The gain is unity power, and L is the pure time delay. The energy balance model under step input is equivalent to a first-order inertial plus pure time delay model. Its unity power gain K, time constant τ, and pure time delay L can be directly identified through step test without the need to measure other physical parameters for derivation.

[0012] Furthermore, based on the equilibrium point power function, the heating and cooling curves are derived, and the identification parameters are obtained, including: From the equilibrium power function, we know that α = 1 / (η) h P max R th ), and the gain K=η h P max R th Therefore, K = 1 / α; Heating direction identification parameters: Based on the parameter model of equivalent verification, the heating response curve T is obtained. h The (t) model is represented as: T h (t)=T h0 +K h ·(u h1 -u h0 )·(1-e -(t-Lh) / τh ), t>L h ; Among them, T h0 U is the initial temperature for the heating experiment. h0 u is the initial power duty cycle for the heating experiment. h1 K represents the final power duty cycle of the heating experiment. h For the heating power gain, τ h Heating time constant, L h Pure lag time for temperature rise; K h =1 / α h ; When calculating the coefficients in the equilibrium point power function using directional methods, the heating direction coefficient α is used. h When using a uniform fitting method for the equilibrium point power function, α is calculated separately according to the direction-specific calculation method. h ; Nonlinear least squares fitting is used to directly optimize τ h L h The temperature response curve T h (t) The sum of squared errors between the model output and the measured curve is minimized; Using a graphical method, the interval from the initial step change to the moment when the temperature first deviates from the initial steady-state value is identified on the heating and cooling curves as L. h And find the temperature rise to (T h1 -T h0 The time t is 0.632 × 10⁻⁶. 0.632 , then τ h =t 0.632 -L h ; Obtain the heating power gain K in the heating direction.h Heating time constant τ h Pure lag time L h ; Cooling direction identification parameters: Similarly, the cooling response curve T c The (t) model is represented as: T c (t)=T c0 +K c ·(u c1 -u c0 )·(1-e -(t-Lc) / τc ), t>L c ; T c0 u is the initial temperature for the cooling experiment. c0 u is the initial power duty cycle for the cooling experiment. c1 K is the final power duty cycle for the cooling experiment. c Cooling power gain, τ c L is the cooling time constant. c This refers to the pure lag time for cooling. Among them, the cooling power gain K c The calculation is the same as the heating power gain, yielding the following results: K c =1 / α c ; The cooling time constant τ is determined using an identification method that aligns with the heating direction. c and cooling pure lag time L c .

[0013] Furthermore, in step 3, the net power deviation is defined as including: In the dynamic simulation, assuming the current simulation step size is Ts, and the current simulation temperature T(k) and current input power u(k) are known; the equilibrium power u at the current temperature is calculated based on the equilibrium point power function from step 1. bal (k); Based on the current input power u(k), and further at the current temperature, the equilibrium power u bal (k) Define the net power deviation Δu(k) as follows: Δu(k)=u(k)-u bal (k); When Δu(k)>0, it indicates that the system is in a state of heating or active heating. When Δu(k)<0, it indicates that the system is in a state of passive heat dissipation or active cooling. When Δu(k)=0, it indicates that the system is in a state of equilibrium, and the heating power just maintains the current temperature.

[0014] Furthermore, in step 3, the identification parameters are smoothly transitioned based on the defined net power deviation to obtain the gain parameters, including: To avoid abrupt changes in model parameters, a smooth transition is achieved for the parameters of the heating and cooling models based on the net power deviation Δu. First, a linear transition interval [-Δu] is defined. max ,Δu max The weighting coefficient λ changes linearly with the net power deviation Δu: When Δu≤-Δu max At that time, λ(Δu) = 0; When -Δu max <Δu<Δu max At that time, λ(Δu) = (Δu + Δu) max ) / (2Δu max ); When Δu≥Δu max At that time, λ(Δu) = 1; Where Δu max The maximum net power deviation threshold is set. Smooth weighting is applied to the time constant τ(Δu) under net power deviation and the pure time delay L(Δu) under net power deviation: τ(Δu)=λ(Δu)τ h +(1-λ(Δu))τ c ; L(Δu) = λ(Δu)L h +(1-λ(Δu))L c ; Under the assumption of ideal linear thermal resistance, the steady-state gain in the direction of system heating and cooling should be equal, i.e., K h =K c =K; In actual temperature control systems, there are thermal resistance direction dependence or cooling nonlinearity factors, which affect the K value identified from step data. h With K c Not completely equal; to balance physical consistency and transient simulation accuracy, the power gain under net power deviation is also smoothly weighted during dynamic simulation: K(Δu) = λ(Δu)K h +(1-λ(Δu))K c .

[0015] Further, in step 3, a first-order linear differential model is established based on the equilibrium point power function, discretized using the first-order backward difference method, and recursively derived based on a set simulation step size to obtain the recursive model, which includes: For any time-varying input u(t), the system is no longer in a single step response mode. Starting from the original energy balance equation, based on the definition of the equilibrium power function in step 1 and the power gain K=η in the equivalent verification in step 2,h P max R th ; The relation u is obtained bal =(TT amb Substituting K / K into the energy balance differential equation, we obtain the first-order linear differential equation for temperature: τ·(dT / dt)+T=K·u(t)+T amb (t); Where τ=R th C, T amb (t) represents the ambient temperature; Discretization is performed using the first-order backward difference method, with a simulation step size of T. s Delay steps d=[L / T s The recurrence relation is expressed as follows: T(k)=a×T(k-1)+b×u(kd)+c(k); in: a=τ / (τ+T s ),b=(K×T s ) / (τ+T s ),c(k)=(T amb (k)×T s ) / (τ+T s ); K, τ, and L are the gain, time constant, and pure time delay calculated in step 3 based on the current net power deviation Δu(k); Ts is the simulation step size; T(k) is the simulated temperature output at step k, and the subscript k is the discretization step index; u(kd) is the input power duty cycle at step kd; d is L / Ts, which is calculated and rounded to the nearest integer, representing the delay step corresponding to pure time lag. T amb c(k) is the ambient temperature at step k. If the ambient temperature is constant, then c(k) is a constant.

[0016] The technical solution adopted in this invention has the following beneficial effects: 1. Reproducing Real-World Dynamics: All parameters (K) of the model in this invention h ,τ h ,L h ,K c ,τ c ,L c The simulation model (including the equilibrium point function) is directly derived from the step response data of the real system, rather than theoretical approximations or empirical assumptions. The established simulation model can more accurately reproduce the asymmetric dynamic behavior of the real temperature control system, providing a reliable simulation method for digital twin of temperature and temperature controller design.

[0017] 2. Simple data acquisition: Only one heating step and one cooling step are needed to obtain all the data required for calculation.

[0018] 3. No physical parameters required: It does not rely on physical quantities that are difficult to measure, such as heat capacity and thermal resistance, making it highly practical for engineering applications.

[0019] 4. High model accuracy: By using two sets of FOPDT models with weighted switching, it can adapt to the dynamic differences in heating and cooling.

[0020] 5. Incorporation of ambient temperature variation: This simulation model supports the inclusion of ambient temperature as a modeling variable. Even in scenarios with significant ambient temperature variations, it can still accurately calculate the inherent steady-state gain of the system, thus expanding the model's applicability.

[0021] 6. Continuous and smooth simulation: Based on linear calculation of net power deviation, the parameters are smoothly transitioned, avoiding curve jumps.

[0022] 7. Strong physical interpretability: This invention derives the equilibrium point power function u through the energy balance equation. bal This makes the net power deviation Δu=uu bal (T) has a clear physical meaning (direction of net heating power). However, the physical meaning of traditional model switching methods is unclear, and this invention effectively overcomes this deficiency.

[0023] 8. Wide applicability: Since heating and cooling asymmetry is a common phenomenon in industrial temperature control systems, this invention can be widely applied to heating-dominated, cooling-dominated, and bidirectional asymmetric systems. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating a simulation modeling method for an asymmetric temperature control system based on energy balance, provided by the present invention. Detailed Implementation

[0025] To make the objectives, technical solutions, and effects of the present invention clearer and more explicit, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0026] This invention provides a simulation modeling method for an asymmetric temperature control system based on energy balance. The simulation modeling method includes: Step 1: Based on the physical system setting of heating power, conduct heating step experiment and cooling step experiment to obtain steady-state temperature and power duty cycle data, and obtain the equilibrium point power function; Step 2: Perform equivalent verification of the equilibrium point power function based on the first-order inertial plus pure time delay model, and obtain the heating curve and cooling curve based on the equilibrium point power function to obtain the identification parameters. Step 3: Define the net power deviation, and based on the defined net power deviation, smooth the transition of the identification parameters to obtain the gain parameters. Establish a first-order linear differential model according to the equilibrium point power function, discretize it using the first-order backward difference method, and recursively deduce based on the set simulation step size to obtain the recursive model. Step 4: Repeat step 3 based on the simulation step size to obtain multiple sets of recursive models, establish the temperature response curve of the temperature control system, and complete the simulation modeling.

[0027] More specifically, in combination Figure 1 The simulation modeling method of the present invention is described in the following process: Step 1: Based on the physical system setting of heating power, conduct heating step experiment and cooling step experiment to obtain steady-state temperature and power duty cycle data, and further obtain the equilibrium point power function; Temperature step test: Set a small heating power in the physical system, for example, u=10%, and record the initial temperature T after steady-state heating. h0 Initial power duty cycle u h0 Initial ambient temperature T amb_h0 Then, the heating power is switched to a higher level, for example, u=100%, and the temperature response curve T of the entire heating process is recorded. h (t), until the system stabilizes again, and record the final steady-state temperature T. h1 Final heating power duty cycle u h1 , steady-state ambient temperature T amb_h1 .

[0028] Cooling step test: Set a relatively large heating power in the physical system, for example, u=100%, and record the initial cooling temperature T after steady-state cooling. c0 Initial power duty cycle for cooling, u c0 Initial ambient temperature T during cooling amb_c0 Then, the heating power is switched to a lower level, for example, u=10%, and the temperature response curve T of the entire cooling process is recorded. c (t), until the system stabilizes again, and record the final cooling steady-state temperature T. c1 Final cooling power duty cycle u c1 Cooling steady-state ambient temperature T amb_c1 .

[0029] Explanation of test data: The two step tests obtained a total of four steady-state points: before heating (T) h0 ,u h0 ,T amb_h0 After heating (T) h1 ,u h1,T amb_h1 ), before cooling (T) c0 ,u c0 ,T amb_c0 (usually u) c0 =u h1 ,T c0 =T h1 ,T amb_c0 =T amb_h1 After cooling (T) c1 ,u c1 ,T amb_c1 (Typically, experimental design will take u) h0 =u c1 This means that the target power duty cycle after cooling is the same as the initial power duty cycle before heating; it should be noted that due to ambient temperature drift or errors in recording the steady-state temperature, the steady-state temperature T after cooling may vary. c1 It may not be equal to the initial steady-state temperature T. h0 Ambient temperature T amb_c1 It may also be related to T amb_h0 (Different). Initial temperature T h0 T c0 The system was already in a stable state before the experiment began, meaning it had been running at a fixed power for a sufficient period of time; the final steady-state temperature T... h1 T c1 If it cannot be strictly achieved, the temperature change rate threshold method (e.g., the absolute value of the temperature change rate over a continuous period of time is less than 0.01°C / s) or the exponential fitting method can be used to determine it.

[0030] Calculate the equilibrium point power function based on the experimental data; Establish the energy balance equation: C·(dT / dt)=η h ·P max ·u-(TT amb ) / R th (1); Where T is the system temperature (°C), t is time (s), C is the total heat capacity of the system (J / °C), and η is the system temperature. h For heating efficiency, P max T is the maximum power of the heater (W), u is the power duty cycle (i.e., the ratio of actual heating power to maximum power, ranging from 0 to 1 or 0% to 100%), and T is the maximum power of the heater (W). amb R represents ambient temperature (°C). th The equivalent thermal resistance (°C / W) between the system and the environment. If the actual system uses a voltage duty cycle u... v Control, based on the principle that heating power is proportional to the square of voltage, can be expressed as u = u v 2 In this invention, a power duty cycle u is used uniformly to maintain the linearity of the model.

[0031] In steady state, (dT / dt=0), substituting into equation (1) yields the equilibrium power u. bal The formula is as follows: u bal =1 / (η h ·P max ·R th )·TT amb / (η h ·P max ·R th (2); Let α = 1 / (η) h ·P max ·R th Then equation (2) can be written as: u bal =αT-αT amb ; The equilibrium point power function model is defined in two cases: (1) Under the condition of constant ambient temperature: Using a univariate linear equation: u bal (T) = αT + γ, at which point theoretically γ = -αT amb However, constraints may not be enforced during fitting.

[0032] (2) Changes in ambient temperature: Using a bivariate linear equation: u bal (T,T amb )=αT+βT amb In this case, β is theoretically equal to -α, but it is not strictly constrained during fitting.

[0033] Method for determining coefficients: Method 1 (calculate by direction and take average): Using the two steady-state points (T) in the heating stage of step S1 h0 ,u h0 ,T amb_h0 ) and (T h1 ,u h1 ,T amb_h1 Substituting these values ​​into the corresponding formal equation, we obtain the coefficient (α) for the heating direction. h γ h ) or (α) h β h Similarly, using the two steady-state points (T) in the cooling stage of step S1... c0 ,u c0 ,T amb_c0 ) and (T c1 ,u c1 ,T amb_c1 Solving for the coefficient (α) in the cooling direction yields the solution. cγ c ) or (α) c β c ).

[0034] For the unary form, take α=(αh+αc) / 2, γ=(γh+γc) / 2; for the binary form, take α=(αh+αc) / 2, β=(βh+βc) / 2.

[0035] Method 2 (Unified Fitting by Least Squares): Substitute all steady-state points into the corresponding function model and solve for the optimal coefficients directly using the least squares method.

[0036] Both methods can yield the power balance point function u. bal (T) or u bal (T,T amb ).

[0037] Step 2: Perform equivalent verification of the equilibrium point power function based on the first-order inertial plus pure time delay model, and derive the heating curve and cooling curve based on the equilibrium point power function to further obtain identification parameters; Consider a step input: at t=0, the heating power jumps from u0 to u1 and remains thereafter. Let the initial steady-state temperature be T0, satisfying u0=ubal(T0). Substituting equation (2) from step S2 into equation (1), and using the initial conditions, the temperature response can be derived as: T(t) = T0 + K·(u1-u0)·(1-e -(t-L) / τ (3); Where τ=R th C is the time constant, K=η h P max R th L represents unity power gain (°C / %), and L represents pure time lag (usually in seconds, consistent with t).

[0038] Expression (3) is completely consistent with the step response form of the classical first-order inertial plus pure time delay (FOPDT) model. Therefore, the energy balance model is equivalent to the FOPDT model under step input, and its parameters K, τ, L can be directly identified through step test without measuring C and η. h P max R th Physical parameters, etc.

[0039] Identify parameters from heating and cooling curves: From equation (2), we know that α = 1 / (η) h P max R th ), and the gain K=η h P max Rth Therefore, K = 1 / α.

[0040] (1) Identification of heating direction parameters: According to formula (3), the heating response curve T h The (t) model is: T h (t)=T h0 +K h ·(u h1 -u h0 )·(1-e -(t-Lh) / τh ), t>L h (4); Among them, T h0 U is the initial temperature for the heating experiment. h0 u is the initial power duty cycle for the heating experiment. h1 K represents the final power duty cycle of the heating experiment. h For the heating power gain, τ h Heating time constant, L h Pure lag time for temperature rise.

[0041] K h =1 / α h ; If directional calculations are used in step S2, then the heating direction coefficient α should be used. h If a uniform fitting method is used, α should be calculated separately according to the method for calculating each direction. h .

[0042] τ h L h Parameter identification method: Nonlinear least squares fitting: Direct optimization of τ h L h This minimizes the sum of squared errors between the model output and the measured curve. Graphical method: Find the interval on the curve from the initial step change moment to the moment when the temperature first deviates from the initial steady-state value as L. h And find the temperature rise to (T h1 -T h0 The time t is 0.632 × 10⁻⁶. 0.632 , then τ h =t 0.632 -L h .

[0043] From this, the heating direction parameter is obtained: heating power gain K. h Heating time constant τ h Pure lag time L h .

[0044] (2) Identification of cooling direction parameters: Similarly, the cooling response curve T c (t) is: T c (t)=T c0 +K c ·(u c1 -u c0 )·(1-e -(t-Lc) / τc ), t>L c (5); T c0 u is the initial temperature for the cooling experiment. c0 u is the initial power duty cycle for the cooling experiment. c1 K is the final power duty cycle for the cooling experiment. c Cooling power gain, τ c L is the cooling time constant. c This refers to the pure lag time for cooling.

[0045] Where the gain K c The calculation is similar: K c =1 / α c ; The cooling time constant τ is determined using an identification method (nonlinear least squares fitting or graphical method) that aligns with the heating direction. c and cooling pure lag time L c .

[0046] Step 3: Define the net power deviation, and based on the defined net power deviation, smooth the transition of the identification parameters to obtain the gain parameters. Establish a first-order linear differential model according to the equilibrium point power function, discretize it using the first-order backward difference method, and recursively deduce based on the set simulation step size to obtain the recursive model. In the dynamic simulation, assuming the current simulation step size is Ts, the current simulation temperature T(k) and the current input power u(k) are known. The equilibrium power u at the current temperature is calculated based on the equilibrium point power function obtained in step S2. bal (k).

[0047] Define net power deviation as: Δu(k)=u(k)-u bal (k); Δu(k)>0 indicates that the system is in a state of heating or active heating, Δu(k)<0 indicates that the system is in a state of passive heat dissipation or active cooling, and when Δu(k)=0, it indicates that the system is in a state of equilibrium, at which point the heating power just maintains the current temperature.

[0048] Calculate model parameters; To avoid abrupt changes in model parameters, a smooth transition of parameters between the heating and cooling models is needed based on the net power deviation Δu. This step presents a simple and practical linear weighting scheme.

[0049] Define a linear transition interval [-Δu] max ,Δu max The weighting coefficient λ varies linearly with Δu: When Δu≤-Δu max At that time, λ(Δu) = 0; When -Δu max <Δu<Δu max At that time, λ(Δu) = (Δu + Δu) max ) / (2Δu max ); When Δu≥Δu max At that time, λ(Δu) = 1; Where Δu max This is the maximum net power deviation threshold that is set.

[0050] Smooth weighting is applied to the time constant τ(Δu) under net power deviation and the pure time delay L(Δu) under net power deviation: τ(Δu)=λ(Δu)τ h +(1-λ(Δu))τ c ; L(Δu) = λ(Δu)L h +(1-λ(Δu))L c ; Under the assumption of ideal linear thermal resistance, the steady-state gain in the direction of system heating and cooling should be equal, i.e., K h =K c =K. Actual temperature control systems may have factors such as thermal resistance direction dependence or cooling nonlinearity, which may affect the K value identified from step data. h With K c Not completely equal. To balance physical consistency and transient simulation accuracy, this invention also applies a smooth weighting to the gain during dynamic simulation: K(Δu) = λ(Δu)K h +(1-λ(Δu))K c .

[0051] Discretized difference equation recursion; For any time-varying input u(t), the system is no longer in a single step response mode, and equation (3) cannot be used directly. Starting from the original energy balance equation (1), based on the definition of the equilibrium point power function in step S2 and the gain K=η in step S3... h P max R th The relation u can be obtained.bal =(TT amb Substituting K into the energy balance differential equation (1), we obtain the first-order linear differential equation for temperature: τ·(dT / dt)+T=K·u(t)+T amb (t)(6); Where τ=R th C, T amb (t) represents the ambient temperature (which can vary over time). This equation is applicable to any time-varying input.

[0052] Discretization is performed using the first-order backward difference method, with a simulation step size of T. s Delay steps d=[L / T s From this, we can obtain the recursive formula: T(k)=a×T(k-1)+b×u(kd)+c(k)(7); in: a=τ / (τ+T s ),b=(K×T s ) / (τ+T s ),c(k)=(T amb (k)×T s ) / (τ+T s ); K, τ, and L are the gain, time constant, and pure time delay calculated in step S6 based on the current net power deviation Δu(k); Ts is the simulation step size; T(k) is the simulated temperature output at step k (the subscript k is the discretization step index, which is different from the gain symbol K (uppercase), so be careful to distinguish it); u(kd) is the input power duty cycle at step kd; d is the rounded L / Ts, which is the delay step corresponding to pure lag; T amb c(k) is the ambient temperature at step k (if the ambient temperature is constant, then c(k) is a constant).

[0053] Step 4: Repeat step 3 based on the simulation step size to obtain multiple sets of recursive models, establish the temperature response curve of the temperature control system, and complete the simulation modeling; repeat step 3 at each simulation step size to obtain the complete temperature response curve.

[0054] Based on the above simulation modeling method, the present invention achieves the following technical effects: 1. The pure energy balance model depends on the physical parameters heat capacity C and thermal resistance R. th Heating efficiency η h Parameters that are difficult to obtain are not practical for engineering applications.

[0055] This invention can collect all the data required for modeling through a single heating step and a single cooling step, without the need to measure specific physical parameters. It is simple to operate and highly practical.

[0056] 2. A single FOPDT model can only describe the characteristics in one direction of heating or cooling, and cannot accurately describe the dynamic process of heating and cooling asymmetric at the same time.

[0057] This invention identifies the parameters (τ) of the heated FOPDT model. h ,L h ,K h ) and cooling FOPDT model parameters (τ) c ,L c ,K c The two sets of parameters are then weighted smoothly based on the net power deviation Δu, so that the model is biased towards the heating model parameters during the heating phase and towards the cooling model parameters during the cooling phase, thereby accurately describing the asymmetric dynamic process.

[0058] 3. The hard switching between the two models lacks a basis, and the temperature jumps near the equilibrium point cause discontinuities in the simulation curves; the ambient temperature is not included, so the influence of ambient temperature changes on the model cannot be eliminated.

[0059] This invention derives the equilibrium point power function u based on energy balance. bal (T) (constant ambient temperature) or u bal (T,T amb (Ambient temperature change), incorporating ambient temperature into the modeling variable; using net power deviation Δu=uu bal As the physical basis for model switching, linear weights are used to smoothly transition the gain, time constant, and pure time delay of the heating and cooling models. Parameters change continuously without temperature jumps caused by hard switching, resulting in smooth and continuous simulation curves that accurately reflect the physical transition characteristics.

[0060] The three features of this invention are interdependent and indispensable: without the double step test, two independent sets of FOPDT parameters cannot be obtained; without the energy balance point power function u... bal Without this, net power deviation cannot be defined; without smooth weighting, curve jumps cannot be avoided.

[0061] Therefore, this invention has the technical advantages of being easy to operate, requiring no physical parameters, providing continuous and smooth simulation, and highly faithfully reproducing asymmetric dynamics, thus providing a reliable and practical simulation method for the design of digital twins and controllers for temperature control systems.

[0062] To further verify the effectiveness of the present invention, the following experiments were conducted: Assume a naturally cooled insulated box equipped only with an electric heater. The heater is controlled by an AC solid-state relay, and the control command is a voltage duty cycle u.v (Value range 0%~100%), actual heating power and u v It is proportional to the square of the power. The model of this invention uses a power duty cycle u = u v 2 As input, to maintain the linearity of the energy balance equation. In the following steps, all power duty cycles u are derived from the voltage duty cycle u. v Calculation yields (u=u) v 2 ).

[0063] 1. Step test and steady-state point extraction: Temperature step: Initial voltage duty cycle u v =10%, corresponding to power duty cycle u h0 =10%, stable temperature T h0 =31.0℃, record the ambient temperature T amb_h0 =20.0℃; Step-to-voltage duty cycle u v =100%, corresponding to the power duty cycle u h1 =100%, final stable temperature T h1 =80.0℃, record the ambient temperature T amb_h1 =21.0℃. The steady-state points (31℃, 10%, 21℃) and (80℃, 100%, 21℃) were obtained.

[0064] Cooling step: Initial voltage duty cycle u v =100%, corresponding to the power duty cycle u c0 =100%, stable temperature T c0 =80.0℃, ambient temperature T amb_c0 =21.0℃; Step-to-voltage duty cycle u v =10%, corresponding to power duty cycle u c1 =10%, final stable temperature T c1 =32℃, ambient temperature T amb_c1 =22℃. The steady-state points were obtained as (80℃, 100%, 21.0℃) and (32℃, 10%, 22.0℃).

[0065] Equilibrium point function coefficients calculation (averaging after fitting in each direction): Let the equilibrium point function be u. bal (T,T amb )=αT+βT amb .

[0066] Substituting the two steady-state points (31,0.1,20) and (80,1.0,21) from the temperature step into the equilibrium function: 0.1 = α h×31+β h ×20; 1.0 = α h ×80+β h ×21; Solving for α, we get: h ≈0.0189, β h ≈-0.0242.

[0067] Substituting the two steady-state points (80, 1.0, 21) and (32, 0.1, 22) from the cooling step into the equilibrium function: 1.0 = α c ×80+β c ×21; 0.1 = α c ×32+β c ×22; Solving for α, we get: c ≈0.0183, β c ≈-0.0221; Take the average value: α=(0.0189+0.0183) / 2=0.0186; β=(-0.0242-0.0221) / 2=-0.02315; The equilibrium point power function is obtained as: u bal (T,T amb = 0.0186T - 0.02315T amb .

[0068] 2. Parameter identification: τ and L are identified from the temperature step curve using a graphical method.

[0069] Suppose that the time constant and lag time are identified from the heating curve as follows: τ h =30s, L h =5s.

[0070] Calculate the gain K h =1.0 / α h =1.0 / 0.0189≈52.91.

[0071] Suppose that the time constant and lag time are identified from the cooling curve as follows: τ c =120s, L c =8s.

[0072] Calculate the gain K c =1.0 / α c=1.0 / 0.0183≈54.64.

[0073] 3. Simulation recursion: Set simulation step size T s =1s, set the transition region Δu max =0.05 (corresponding to a 5% power change). The initial state of the simulation system is set as follows: initial temperature T(0) = 30℃, initial power u(0) = 0%. The power output is set as follows (in actual control, the voltage duty cycle is used; the power duty cycle needs to be square rooted and converted to voltage duty cycle): 0–100s: Power increases uniformly from 0% to 50% (corresponding to the voltage duty cycle increasing from 0% to approximately 70.71% according to the square root curve of power); 100–200s: Power maintained at 50% (voltage duty cycle approximately 70.71%); 200–300s: Power decreases uniformly from 50% to 0% (voltage duty cycle decreases from 70.71% to 0% according to the square root curve of power); 300–400s: Power remains at 0% (voltage duty cycle 0%); 400–500s: Power increases uniformly from 0% to 100% (voltage duty cycle increases from 0% to 100% according to the square root curve of power).

[0074] 500-600s: Power maintained at 100% (voltage duty cycle 100%).

[0075] At each simulation step, calculate the current net power deviation Δu(k) = u(k) - u bal (T(k)), calculate the coefficient λ according to the linear weight formula, and obtain the gain K, time constant τ and pure lag L of the current step size by weighting. Substitute them into the recursive formula (7) to calculate the temperature at the next moment.

[0076] The simulation started at t=0 and was repeated second by second until t=600, yielding the temperature response curve. The simulation curve exhibited typical asymmetric characteristics: a rapid increase in temperature during the heating phase, a slow decrease during the cooling phase, and a smooth transition without jumps during power abrupt changes, consistent with the expected dynamic behavior, thus verifying the effectiveness of the model established in this invention.

Claims

1. A simulation modeling method for an asymmetric temperature control system based on energy balance, characterized in that, Simulation modeling methods include: Step 1: Based on the physical system setting of heating power, conduct heating step experiment and cooling step experiment to obtain steady-state temperature and power duty cycle data, and obtain the equilibrium point power function; Step 2: Perform equivalent verification of the equilibrium point power function based on the first-order inertial plus pure time delay model, and obtain the heating curve and cooling curve based on the equilibrium point power function to obtain the identification parameters. Step 3: Define the net power deviation, and based on the defined net power deviation, smooth the transition of the identification parameters to obtain the gain parameters. Establish a first-order linear differential model according to the equilibrium point power function, discretize it using the first-order backward difference method, and recursively deduce based on the set simulation step size to obtain the recursive model. Step 4: Repeat step 3 based on the simulation step size to obtain multiple sets of recursive models, establish the temperature response curve of the temperature control system, and complete the simulation modeling.

2. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, In step 1, the temperature step experiment includes: In the physical system, a low heating power u = 10% is set, and the initial temperature T after steady-state heating is recorded. h0 Initial power duty cycle u h0 Initial ambient temperature T amb_h0 Then, the heating power is switched to a high heating power u=100%, and the temperature response curve T of the entire heating process is recorded. h (t), until the system stabilizes again, and record the final steady-state temperature T. h1 Final heating power duty cycle u h1 , steady-state ambient temperature T amb_h1 ; The cooling step test includes: In the physical system, a high heating power u=100% is set, and the initial temperature T after cooling down after steady state is recorded. c0 Initial power duty cycle for cooling u c0 Initial ambient temperature T during cooling amb_c0 Then, the heating power is switched to a low heating power u=10%, and the temperature response curve T of the entire cooling process is recorded. c (t), until the system stabilizes again, and record the final cooling steady-state temperature T. c1 Final cooling power duty cycle u c1 Cooling steady-state ambient temperature T amb_c1 .

3. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, In step 1, based on the steady-state temperature and power duty cycle data, the equilibrium point power function is obtained, and the energy balance equation is first established: C·(dT / dt)=η h ·P max ·u-(T-T amb ) / R th ; Where T is the system temperature, t is time, C is the total heat capacity of the system, and η is the system temperature. h For heating efficiency, P max T represents the heater's maximum power, u is the power duty cycle (the ratio of actual heating power to maximum power), and its value ranges from 0 to 1 or from 0% to 100%. amb R represents the ambient temperature. th The equivalent thermal resistance between the system and the environment; In steady state, dT / dt = 0. Substituting this into the energy balance equation, we obtain the equilibrium power u. bal The formula is as follows: you bal =1 / (n h ·P max ·R th )·TT amb / (or h ·P max ·R th ); Let α = 1 / (η) h ·P max ·R th If the equilibrium point power u is... bal Represented as: u bal =αT-αT amb ; When the ambient temperature is constant, the equilibrium point power function is expressed as a univariate linear equation as follows: the bal (T)=αT+γ; When the ambient temperature changes, the equilibrium power function is expressed as a two-variable linear equation: u bal (T,T amb )=αT+βT amb 。 4. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 3, characterized in that, The coefficients in the equilibrium point power function are obtained using either a directional calculation method or the least squares method. The directional calculation includes: Using the two steady-state points from the temperature increase step experiment in step 1, substitute them into the corresponding formal equation to obtain the coefficient (α) for the heating direction. h γ h ) or (α) h β h Similarly, the coefficient (α) of the cooling direction can be obtained by solving the two steady-state points of the cooling step test in step 1. c γ c ) or (α) c β c ); For the unary form, take α = (α h +α c ) / 2, γ=(γ h +γ c ) / 2; For the binary form, take α=(α h +α c ) / 2, β=(β h +β c ) / 2; Least squares methods include: Substitute all steady-state points into the corresponding function model and solve for the optimal coefficients directly using the least squares method.

5. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, Step 2, the equivalent verification of the equilibrium power function based on the first-order inertial plus pure time delay model includes: Considering a step input, at t=0, the heating power jumps from the initial power u0 to the final power u1 and remains thereafter; let the initial steady-state temperature be T0, satisfying u0=u bal (T0); Based on the equilibrium point power model from step 1, and using the initial conditions, the temperature response is derived as follows: T(t)=T0+K·(u1-u0)·(1-e -(t-L) / τ ); Where τ=R th C is the time constant, K=η h P max R th The gain is unity power, and L is the pure time delay. The energy balance model under step input is equivalent to a first-order inertial plus pure time delay model. Its unity power gain K, time constant τ, and pure time delay L can be directly identified through step test without the need to measure other physical parameters for derivation.

6. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, Based on the equilibrium point power function, the heating and cooling curves are derived, and the identification parameters include: From the equilibrium power function, we know that α = 1 / (η) h P max R th ), and the gain K=η h P max R th Therefore, K = 1 / α; Heating direction identification parameters: Based on the parameter model of equivalent verification, the heating response curve T is obtained. h The (t) model is represented as: T h (t)=T h0 +K h ·(u h1 -u h0 )·(1-e -(t-Lh) / τh ),t>L h ; Among them, T h0 u is the initial temperature for the heating experiment. h0 u is the initial power duty cycle for the heating experiment. h1 K represents the final power duty cycle of the heating experiment. h For the power gain due to heating, τ h Heating time constant, L h Pure lag time for temperature rise; K h =1 / α h ； When calculating the coefficients in the equilibrium point power function using directional methods, the heating direction coefficient α is used. h When using a uniform fitting method for the equilibrium point power function, α is calculated separately according to the direction-specific calculation method. h ; Nonlinear least squares fitting is used to directly optimize τ h L h The temperature response curve T h (t) The sum of squared errors between the model output and the measured curve is minimized; Using a graphical method, the interval from the initial step change to the moment when the temperature first deviates from the initial steady-state value is identified on the heating and cooling curves as L. h And find the temperature rise to (T h1 -T h0 The time t is 0.632 × 10⁻⁶. 0.632 , then τ h =t 0.632 -L h ; Obtain the heating power gain K in the heating direction. h Heating time constant τ h Pure lag time L h ; Cooling direction identification parameters: Similarly, the cooling response curve T c The (t) model is represented as: T c (t)=T c0 +K c ·(u c1 -u c0 )·(1-e -(t-Lc) / τc ),t>L c ; T c0 u is the initial temperature for the cooling experiment. c0 u is the initial power duty cycle for the cooling experiment. c1 K is the final power duty cycle for the cooling experiment. c Cooling power gain, τ c L is the cooling time constant. c This refers to the pure lag time for cooling. Among them, the cooling power gain K c The calculation is the same as the heating power gain, yielding the following results: K c =1 / α c ； The cooling time constant τ is determined using an identification method that aligns with the direction of heating. c and cooling pure lag time L c .

7. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, In step 3, the net power deviation is defined as including: In the dynamic simulation, assuming the current simulation step size is Ts, and the current simulation temperature T(k) and current input power u(k) are known; the equilibrium power u at the current temperature is calculated based on the equilibrium point power function from step 1. bal (k); Based on the current input power u(k), and further at the current temperature, the equilibrium power u bal (k) Define the net power deviation Δu(k) as follows: Δu(k)=u(k)-u bal (k); When Δu(k)>0, it indicates that the system is in a state of heating or active heating. When Δu(k)<0, it indicates that the system is in a state of passive heat dissipation or active cooling. When Δu(k)=0, it indicates that the system is in a state of equilibrium, and the heating power just maintains the current temperature.

8. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, In step 3, the identification parameters are smoothly transitioned based on the defined net power deviation to obtain the gain parameters, including: To avoid abrupt changes in model parameters, a smooth transition is achieved for the parameters of the heating and cooling models based on the net power deviation Δu. First, a linear transition interval [-Δu] is defined. max ,Δu max The weighting coefficient λ changes linearly with the net power deviation Δu: When Δu ≤ -Δu max then λ(Δu) = 0; When -Δu max <Δu<Δu max When, λ(Δu)=(Δu + Δu max ) / (2Δu max ); When Δu≥Δu max At that time, λ(Δu) = 1; Where Δu max The maximum net power deviation threshold is set. Smooth weighting is applied to the time constant τ(Δu) under net power deviation and the pure time delay L(Δu) under net power deviation: τ(Δu)=λ(Δu)τ h +(1-λ(Δu))τ c ; L(Δu)=λ(Δu)L h +(1-λ(Δu))L c ; Under the assumption of ideal linear thermal resistance, the steady-state gain in the direction of system heating and cooling should be equal, i.e., K h =K c =K; In actual temperature control systems, factors such as thermal resistance direction dependence or cooling nonlinearity may affect the K value identified from step data. h With K c Not completely equal; to balance physical consistency and transient simulation accuracy, the power gain under net power deviation is also smoothly weighted during dynamic simulation: K(Δu)=λ(Δu)K h +(1-λ(Δu))K c 。 9. The simulation modeling method for an asymmetric temperature control system based on energy balance according to claim 1, characterized in that, In step 3, a first-order linear differential model is established based on the equilibrium point power function, discretized using the first-order backward difference method, and recursively derived based on a set simulation step size to obtain the recursive model, which includes: For any time-varying input u(t), the system is no longer in a single step response mode. Starting from the original energy balance equation, based on the definition of the equilibrium point power function in step 1 and the power gain K=η in the equivalent verification in step 2, h P max R th ; The relation u is obtained bal =(TT amb Substituting K / K into the energy balance differential equation, we obtain the first-order linear differential equation for temperature: τ·(dT / dt)+T=K·u(t)+T amb (t); Where τ=R th C, T amb (t) represents the ambient temperature; Discretization is performed using the first-order backward difference method, with a simulation step size of T. s Delay steps d=[L / T s The recurrence relation is expressed as follows: T(k)=a×T(k-1)+b×u(kd)+c(k); in: a=τ / (τ+T s ),b=(K×T s ) / (τ+T s ),c(k)=(T amb (k)×T s ) / (τ+T s ); K, τ, and L are the gain, time constant, and pure time delay calculated in step 3 based on the current net power deviation Δu(k); Ts is the simulation step size; T(k) is the simulated temperature output at step k, and the subscript k is the discretization step index; u(kd) is the input power duty cycle at step kd; d is L / Ts, which is calculated and rounded to the nearest integer, representing the delay step corresponding to pure time lag. T amb c(k) is the ambient temperature at step k. When the ambient temperature is constant, c(k) is a constant.

Images