CN113158546B - Fault prediction method of air conditioning unit equipment based on HSMM - Google Patents

Fault prediction method of air conditioning unit equipment based on HSMM Download PDF

Info

Publication number
CN113158546B
CN113158546B CN202110262737.0A CN202110262737A CN113158546B CN 113158546 B CN113158546 B CN 113158546B CN 202110262737 A CN202110262737 A CN 202110262737A CN 113158546 B CN113158546 B CN 113158546B
Authority
CN
China
Prior art keywords
state
time
cooling coil
fan
fault
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110262737.0A
Other languages
Chinese (zh)
Other versions
CN113158546A (en
Inventor
严颖
蔡骏
张菀
李涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN202110262737.0A priority Critical patent/CN113158546B/en
Publication of CN113158546A publication Critical patent/CN113158546A/en
Application granted granted Critical
Publication of CN113158546B publication Critical patent/CN113158546B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/29—Graphical models, e.g. Bayesian networks
    • G06F18/295—Markov models or related models, e.g. semi-Markov models; Markov random fields; Networks embedding Markov models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/084—Backpropagation, e.g. using gradient descent
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02B—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO BUILDINGS, e.g. HOUSING, HOUSE APPLIANCES OR RELATED END-USER APPLICATIONS
    • Y02B30/00—Energy efficient heating, ventilation or air conditioning [HVAC]
    • Y02B30/70—Efficient control or regulation technologies, e.g. for control of refrigerant flow, motor or heating

Abstract

The invention relates to a fault prediction method of air conditioning unit equipment based on HSMM, firstly adopting methods such as neural network and filtering to extract fault characteristics, and combining a principal component analysis method to obtain a principal fault characteristic matrix; then, establishing initial semi-hidden Markov models corresponding to faults of the air conditioning unit respectively, and training by using a Gibbs sampling method to obtain the semi-hidden Markov models; then, a scaling Viterbi algorithm is applied to obtain state estimation of a semi-hidden Markov model corresponding to each device of the air conditioning unit; then, a backward recursion algorithm is applied to synthesize the semi-hidden Markov models of different faults and the semi-hidden Markov models of the air conditioning unit and equipment thereof, so that the calculated amount is reduced; finally, based on the acquisition of the current state estimation value of the semi-hidden Markov model corresponding to the air conditioning unit and the equipment thereof, further adopting a backward recursion algorithm to calculate and acquire the residual service life of the air conditioning unit and the equipment thereof; therefore, accurate prediction can be realized for the faults of the air conditioning unit equipment, and the actual working efficiency of the air conditioning unit is ensured.

Description

Fault prediction method of air conditioning unit equipment based on HSMM
Technical Field
The invention relates to a fault prediction method of air conditioning unit equipment based on an HSMM, and belongs to the technical field of air conditioning unit fault diagnosis.
Background
About 50% of the energy in commercial buildings is consumed by hvac systems. In hvac systems, faults may cause significant energy consumption and may not meet the needs of people for comfort. As an important subsystem of a hvac system, an air conditioning unit (air handling unit, AHU) is used to condition indoor air to a suitable temperature. The system consists of a mixing box, a filter, a coil pipe, a fan, an air supply pipeline and the like. Failure of the air conditioning unit can affect heat exchange between chilled water and air or air flow, so that energy waste and comfort of people cannot be met. The fault prediction may predict the time and location at which the fault may occur by estimating the remaining useful life (remaining useful life, RUL) of the device before the fault occurs. The fault which causes great harm and suddenly occurs in the future can be avoided by timely maintaining or replacing potential fault equipment, the service life of the equipment is prolonged, and the operation and maintenance cost is reduced. Therefore, it is very critical to accurately and timely perform fault prediction, and the current research on fault prediction of air conditioning units is very little.
Disclosure of Invention
The invention aims to solve the technical problem of providing a fault prediction method of air conditioning unit equipment based on an HSMM, which adopts a brand new design architecture, and adopts a semi-hidden Markov model and a backward recursion algorithm, so that high-precision fault diagnosis prediction can be realized for the air conditioning unit, and the stability of actual work is ensured.
The invention adopts the following technical scheme for solving the technical problems: the invention designs a fault prediction method of air conditioning unit equipment based on HSMM, which comprises the following steps:
step A, based on the time sequence changes of the related variables of the faults of the air conditioning unit, constructing time sequence changes of the fault characteristics corresponding to the cooling coil and the fan on the air conditioning unit, and further obtaining a fault characteristic matrix X corresponding to the cooling coil and the fan cc 、X sf And further aims at the fault characteristic matrix X cc 、X sf Respectively executing fault characteristic dimension reduction operation to obtain main fault characteristic matrixes O corresponding to the cooling coils and the fans respectively cc 、O sf Then enter step B;
step B, according to the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf Constructing initial semi-hidden Markov models corresponding to different faults of the cooling coil respectively, constructing initial semi-hidden Markov models corresponding to different faults of the fan respectively, and based on a main fault characteristic matrix O cc and Osf Training the preset proportion training data of the cooling coil by using a Gibbs sampling method aiming at each initial semi-hidden Markov model respectively to obtain semi-hidden Markov models corresponding to different faults of the cooling coil and semi-hidden Markov models corresponding to different faults of the fan respectively, and then entering the step C;
step C, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf The preset proportion test data in the cooling coil is obtained by adopting a scaling Viterbi algorithm, and the semi-hidden Markov model corresponding to different faults of the cooling coil and different faults of the fan respectively and the state corresponding to the time 1 to the time K are obtainedEstimating states, and then entering a step D;
step D, integrating the semi-hidden Markov models corresponding to the faults of the equipment respectively by applying a backward recursion algorithm to the cooling coil and the fan to obtain an air conditioning unit, the semi-hidden Markov models corresponding to the cooling coil and the fan respectively, further obtaining the current state estimation values of the air conditioning unit, the cooling coil and the fan corresponding to the semi-hidden Markov models, and then entering the step E;
and E, calculating to obtain the residual service lives of the air conditioning unit, the cooling coil and the fan by adopting a backward recursion algorithm according to the air conditioning unit, the semi-hidden Markov model corresponding to the cooling coil and the fan and the current state estimated value of each device.
As a preferred embodiment of the present invention, the step a includes the following steps A1 to A5:
step A1. Collecting the flow of refrigerating water of an air conditioning unit
Figure GDA0004143077470000021
Air flow rate->
Figure GDA0004143077470000022
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA0004143077470000023
And->
Figure GDA0004143077470000024
Respectively used as fault characteristic time sequence changes on the air conditioning unit;
at the same time, the air supply temperature T of the air conditioning unit is collected a,sup And a preset value T thereof sup,spt The difference corresponds to the time sequence change delta T from time 1 to time K sup,spt (1)、…、ΔT sup,spt (K) The method comprises the steps of taking the fault characteristic time sequence change as a fault characteristic time sequence change on an air conditioning unit, and then entering a step A2;
a2, constructing the refrigerating water flow corresponding to the training cooling coil
Figure GDA0004143077470000025
Flow of mixed wind->
Figure GDA0004143077470000026
Temperature T of mixed air a,mix Delivery temperature T with chilled water chw,sup The difference is used as input, and the estimated value of the heat exchange quantity of the cooling coil is +.>
Figure GDA0004143077470000027
Reverse-transfer neural network for output and application +.>
Figure GDA0004143077470000028
T a,mix -T chw,sup Corresponding to time sequence changes from time 1 to time K respectively, obtaining a heat exchange quantity estimated value +.>
Figure GDA0004143077470000029
Corresponding to the time sequence change from time 1 to time K;
then according to the specific heat capacity c of the chilled water pw And the chilled water return temperature T chw,rn The following formula applies:
Figure GDA0004143077470000031
obtaining the heat exchange quantity Q cc Corresponding to the time sequence change from time 1 to time K, and then
Figure GDA0004143077470000032
Obtaining DeltaQ cc Time sequence change DeltaQ corresponding to time 1 to time K cc (1)、…、ΔQ cc (K) As the time sequence change of the fault characteristic on the air conditioning unit, entering into the step A3;
A3, obtaining an estimated value of the inner diameter of the upper pipe wall of the cooling coil pipe on the air conditioning unit by applying an ash box model and a particle filtering method
Figure GDA0004143077470000033
Estimation of the area of the ribCount->
Figure GDA0004143077470000034
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA0004143077470000035
And->
Figure GDA0004143077470000036
Obtaining an efficiency estimate of the fan +.>
Figure GDA0004143077470000037
Corresponding to the time sequence change from time 1 to time K
Figure GDA0004143077470000038
Then enter step A4;
a4, obtaining fault characteristic matrixes X corresponding to the cooling coil and the fan respectively cc 、X sf The method comprises the following steps:
Figure GDA0004143077470000039
Figure GDA00041430774700000310
then enter step A5;
step A5. For the failure feature matrix X cc 、X sf Respectively executing PCA principal component analysis method, and aiming at fault feature dimension reduction operation, obtaining a principal fault feature matrix O corresponding to the cooling coil and the fan respectively cc 、O sf 。
As a preferred embodiment of the present invention, the step A3 includes the following steps A3-1 to A3-6:
step A3-1, constructing a cooling coil corresponding to the inner diameter d of the tube wall tube,in Area A of the rib fin The ash bin model of (2) is as follows:
Figure GDA00041430774700000311
wherein ,Ea,mix ＝(1.006T a,mix )+W a,mix (1.84T a,mix +2501),
E a,sup ＝(1.006T a,sup )+W a,sup (1.84T a,sup +2501),
LMTD＝(ΔT sup -ΔT rn )/(lnΔT sup -lnΔT rn )，
τ＝(A fin +A tube,out )/A tube,in ，
α a,e ＝α a [1-A fin (1-η fin )/(A tube,out +A fin )]，
Figure GDA0004143077470000041
wherein ,Wa,mix and Wa,sup Respectively showing the humidity of the mixed air of the air conditioning unit, the humidity of the cooled air supply and E a,mix and Ea,sup Respectively representing the enthalpy value of mixed air of an air conditioning unit and the enthalpy value delta of cooled air supply tube Represents the thickness of the cooling coil pipe, lambda tube Representing the coefficient of thermal conductivity of the cooling coil pipe, R f Indicating the thermal resistance of the cooling coil due to dust, A tube,in Indicating the inner surface area of the cooling coil pipe, A tube,out Represents the outer surface area, deltaT, of the cooling coil pipe sup ＝T a,sup -T chw,sup ,ΔT rn ＝T a,mix -T chw,rn ，T a,sup T represents the air supply temperature chw,rn Represents the return water temperature of chilled water in the cooling coil, alpha a Representing the convective heat transfer coefficient, eta of the rib fin The fin efficiency is represented, the physical coefficient is represented by A,
Figure GDA0004143077470000042
represents the cooling coil cooling water flow rate, and ψ represents the hot fluid densityA degree; then enter step A3-2;
step A3-2. Inner diameter d of pipe wall tube,in Area A of the rib fin As states in the filter frame, a particle filter frame corresponding to the cooling coil is constructed as follows:
cooling coil state equation: x is x cc (t+1)＝x cc (t)+v cc (t)
wherein ,xcc (t+1) represents the state of the cooling coil at the time t+1, x cc (t) represents the state of the cooling coil corresponding to the time t, and the state of the cooling coil comprises the inner diameter d of the tube wall tube,in Area A of the rib fin ，v cc (t) represents the process noise of the cooling coil corresponding to time t, the process noise comprising the process noise of the inner diameter of the tube wall and the process noise of the fin area;
cooling coil measurement equation: z cc (t)＝h(x cc (t))+w cc (t)
wherein ,
Figure GDA0004143077470000043
Figure GDA0004143077470000044
wherein ,wcc (t) representing the measured noise of the cooling coil corresponding to the time t, and then entering the step A3-3;
step A3-3, based on the cooling coil state equation and the cooling coil measurement equation, using a particle filtering method to represent the posterior probability of a random event by a group of weighted random samples, namely particles, and obtaining the estimated value of the inner diameter of the upper tube wall of the cooling coil on the air conditioning unit from the sensor readings containing noise in the formula
Figure GDA0004143077470000051
And rib area estimation +.>
Figure GDA0004143077470000052
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA0004143077470000053
And->
Figure GDA0004143077470000054
Then enter step A3-4;
step A3-4, constructing a fan to include the efficiency e sf The ash bin model of (2) is as follows:
Figure GDA0004143077470000055
wherein ,Qsf Representing the power of the fan; ρ air Representing the density of air;
Figure GDA0004143077470000056
a preset value representing the air supply quantity; ΔP sf,des Representing the pressure drop of the fan, c 1 、c 2 、c 3 、c 4 、c 5 Coefficients of each term in the polynomial, which can be found in a manual of fans, or obtained based on data training using a nonlinear least squares method; then enter step A3-5;
step A3-5, constructing a particle filter frame corresponding to the fan as follows:
fan state equation: e, e sf (t+1)＝e sf (t)+v sf (t)
wherein ,esf (t+1) represents the efficiency of the fan at time t+1, v sf (t) represents process noise of the fan corresponding to the time t;
the fan measurement equation is as follows:
Figure GDA0004143077470000057
wherein ,
Figure GDA0004143077470000058
w sf (t) represents the measurement of the fan at the time tNoise, Q sf (t) represents the power of the fan at the time corresponding to t+1; then enter step A3-6;
step A3-6, based on the fan state equation and the fan measurement equation, obtaining an efficiency estimation value of the fan corresponding to the time t+1 from the noise measurement value and a weighted average of the fan efficiency estimation corresponding to the time t by adopting a particle filtering method
Figure GDA0004143077470000059
Further obtain the efficiency estimation of the fan>
Figure GDA00041430774700000510
Time sequence variation corresponding to time 1 to time K>
Figure GDA00041430774700000511
As a preferred embodiment of the present invention, the step B includes the following steps B1 to B3;
step B1 based on the main fault characteristic matrix O corresponding to the cooling coil cc As observation, determining initial semi-hidden Markov models corresponding to different faults of the cooling coil, wherein each initial semi-hidden Markov model has four states and corresponds to normal, initial fault, middle fault and late fault respectively;
at the same time, based on the main fault characteristic matrix O corresponding to the fan sf As observation, determining initial semi-hidden Markov models corresponding to different faults of the fan respectively, wherein each initial semi-hidden Markov model has four states and corresponds to normal, initial fault, middle fault and late fault respectively; then enter step B2;
b2, respectively constructing a state transition matrix P= { a for each initial semi-hidden Markov model corresponding to different faults of the cooling coil and different faults of the fan based on each semi-hidden Markov model having four states ij },i＝1、…、4,j＝1、…、4，{a ij -representing the probability of the device transitioning from the ith state to the jth state; based on the observed value O can To construct likelihood functions b for observations of different states i (O (s)), i=1,..4, indicating the probability that the observed value O(s) belongs to the state i, then proceeding to step B3;
b3, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf And (3) training the initial semi-hidden Markov models by using a Gibbs sampling method according to the preset proportion training data of the cooling coil, combining the initial semi-hidden Markov models, the state transition matrix and the likelihood functions of the observed values, which correspond to different faults respectively, so as to obtain the semi-hidden Markov models corresponding to the different faults of the cooling coil and the different faults of the fan respectively, and then entering the step C.
As a preferred embodiment of the present invention, the different malfunctions of the cooling coil include the occurrence of a deposit f in the cooling coil pipe tube,cc Dust backlog f on cooling coil rib fin,cc The chilled water valve of the cooling coil is blocked in the full-open state f vlvso,cc The chilled water valve of the cooling coil is blocked in the full-closed state f vlvsc,cc The method comprises the steps of carrying out a first treatment on the surface of the The different faults of the fan comprise the reduction f of the efficiency of the fan ef,sf The fan being locked at a certain speed f fs,sf The fan is thoroughly damaged f cf,sf 。
In the step C, for each semi-hidden markov model corresponding to different faults of the cooling coil and different faults of the fan, according to the following steps C1 to C8, state estimation from time 1 to time K corresponding to the semi-hidden markov model is obtained, that is, state estimation from time 1 to time K corresponding to the semi-hidden markov model corresponding to different faults of the cooling coil and different faults of the fan is obtained;
Step c1, when the initialization time t=1, the hidden state is the probability maximum value in all possible state transition paths of i e {1, 2, …, N };
step C2. forward procedure: the following formula is adopted:
Figure GDA0004143077470000061
calculating the probability delta that the observed value corresponds to the state i of the equipment from the time t-d+1 to the time t i (t) wherein,
Figure GDA0004143077470000071
p i (d) Representing the probability that the device is in state i for d time units; b i (O (s)) represents the probability that the observed value O(s) belongs to the state i;
step C3. is according to the following formula:
Figure GDA0004143077470000072
obtaining
Figure GDA0004143077470000073
wherein />
Figure GDA0004143077470000074
And is based on
Figure GDA0004143077470000075
Calculating a scaling factor->
Figure GDA0004143077470000076
wherein ,bi (O (t)) represents an observation likelihood function corresponding to a time t in the semi-hidden Markov model; when k (t) approaches ≡, the scaling factor k (t) is set to 1 +.>
Figure GDA0004143077470000077
If k (t) is not approaching ≡, reserving the value corresponding to k (t);
step C4, adjusting various parameters in the semi-hidden Markov model based on a scaling factor k (t);
step C5. obtains the maximized probability
Figure GDA0004143077470000078
Duration d of (2), where d=1, …, M i ，M i Is the duration of state iLimiting;
step C6. obtains maximized delta i (t)a ij Corresponding state i, i.e
Figure GDA0004143077470000079
wherein ,aij Representing a probability that the device transitions from the ith state to the jth state;
step C7. backward procedure: calculating the probability of all states at time K and determining the most probable state delta i (K)；
Step C8. obtains the state estimation of the failure semi-hidden Markov model from time 1 to time K by backtracking the recorded d and i
Figure GDA00041430774700000710
In the step D, the following steps D1 to D6 are executed for the cooling coil, the fan and the air conditioning unit, respectively, to obtain the semi-hidden markov model and the current state estimation value corresponding to the air conditioning unit and the cooling coil and the fan thereof, respectively;
d1, generating N according to state transition matrix of semi-hidden Markov model corresponding to different faults of equipment by using Monte Carlo simulation smp A group state sequence, according to the state sequences, randomly selecting a value as the value of the fault parameter in the interval of the fault parameter corresponding to the discrete state, substituting the value of the fault parameter into the ash box model of the cooling coil and the ash box model of the fan, calculating the corresponding equipment capacity, wherein the equipment capacity is the state of the equipment, and further obtaining the fault state s f And device state s dev Sampling points of the relationship between the two; associating each device with a fault state s f Together fault co-operating state s js To indicate that the characterization of the fault joint state s can be further obtained js And device state s dev Sampling points of the relationship between the two;
step D2. is based on the semi-hidden Markov model corresponding to the different faults of the equipment and the state estimation from time 1 to time KThe sampling points obtained in the step D1 are respectively extracted and obtained for each device to obtain two mapping matrixes for reflecting the fault joint state s js And device state s dev A relationship between;
the fault joint state s is embodied js And device state s dev Two mapping matrices of the relationship between each other, the first mapping matrix being represented by a given s dev Time s js Probability of each state, the matrix is defined as M js|dev ＝{m ij,js|dev}, wherein ,
m ij,js|dev ＝Pr(s js ＝j|s dev ＝i)＝n ij,dev /n i,dev ；
in the formula,sdev Is the status of the device; s is(s) js Is the joint state of the fault to which the device corresponds; n is n ij,dev Is corresponding to s dev =i and s js Number of sampling points of =j, n i,dev Is corresponding to state s dev Number of points=i;
corresponding second mapping matrix M dev|js ＝{m ji,dev|js The expression of a given s js Time s dev Probability of each state, where m ji,dev|js ＝Pr(s dev ＝i|s js ＝j)；
Step D3, integrating parameters of the semi-hidden Markov models corresponding to the faults respectively according to the two mapping matrixes of each device to obtain the semi-hidden Markov model corresponding to the device;
step D4., applying a backward recursion algorithm according to the distribution obeyed by the state durations corresponding to different faults respectively, and obtaining the duration of each state of each device under each fault respectively by calculating the evolution time among different joint states;
In step D4, the joint state j belongs to the state k of the device, and when the digital sum of the joint state j is D max,k ：
Figure GDA0004143077470000081
wherein dmax,k Is the maximum value of the sum of the digits of the joint state corresponding to the equipment state k;
Figure GDA0004143077470000082
is a set of joint states in which the sum of digits of the joint states is a maximum;
when the digital sum of the joint state j is d max,k -1：
Figure GDA0004143077470000091
…
When the digital sum of the joint state j is d min,k ：
Figure GDA0004143077470000092
wherein Dj,js Is the duration of the joint state j;
Figure GDA0004143077470000093
from the joint state j to the joint state +.>
Figure GDA0004143077470000094
Is a time of evolution of (2); when->
Figure GDA0004143077470000095
When (I)>
Figure GDA0004143077470000096
The duration of time that is considered to be the device state k; will->
Figure GDA0004143077470000097
The number of joint states of (a) is defined as n min,k The method comprises the steps of carrying out a first treatment on the surface of the If n min,k =1, the duration of device state k is +.>
Figure GDA0004143077470000098
If n min,k >1, then the average value of the durations corresponding to the different joint states j>
Figure GDA0004143077470000099
Can be considered as the duration of state k; />
Step D5. according to embodying the Fault Joint State s js And device state s dev The two mapping matrixes of the relation between the two mapping matrixes, the duration time of each state of the equipment under each fault and the transition matrix of the joint state respectively generate a state transition matrix of a semi-hidden Markov model corresponding to the equipment, and the state transition matrix is obtained according to the following formula for the equipment;
Figure GDA00041430774700000910
wherein ,sjs Is the joint state of the devices; p(s) dev (t+1)＝j|s js (t+1) =n) and P(s) js (t)＝m|s dev (t) =i) from the mapping matrix;
step D6., obtaining the current state estimation value of the semi-hidden Markov model corresponding to the equipment according to the states of different faults of the semi-hidden Markov model corresponding to the equipment;
wherein, based on the three gradual faults of the device being independent, pr(s) js (t)＝i):
Figure GDA00041430774700000911
wherein ,fl ,l＝1,...,N dev Is N of the device dev A fault; joint state s js (t) =i corresponds to N dev Failure f l ,l＝1,...,N dev In state i l ,l＝1,...,N dev In obtaining the joint state s js After the probability of each state, the state of the device at time t is estimated as follows
Figure GDA00041430774700000912
Based on the mapping matrix m ij,dev|js To estimate the number of times the sample is to be processed,
Figure GDA0004143077470000101
wherein Njs Is the number of joint states, and then obtains the current state estimation value of the equipment based on the states of different faults.
As a preferred technical scheme of the invention: in the step E, according to the air conditioning unit, the semi-hidden Markov model corresponding to the cooling coil and the fan of the air conditioning unit and the corresponding current state estimation value, a backward recursion algorithm is adopted, and the remaining service lives of the air conditioning unit, the cooling coil and the fan of the air conditioning unit are obtained through calculation in the following manner;
the air conditioning unit and the semi-hidden Markov model corresponding to the air conditioning unit and the equipment respectively have four states, and the method is implemented as follows, wherein state 0 corresponds to a normal state, and states 1, 2 and 3 correspond to three fault levels respectively;
State 3 RUL 3 ＝0,
State 2 RUL 2 ＝a 22 ·(D 2 -t stay,2 ),
State 1 RUL 1 ＝a 11 ·(D 1 -t stay,1 +RUL 2 )+a 12 ·RUL 2 ,
State 0 RUL 0 ＝a 00 ·(D 0 -t stay,0 +RUL 1 )+a 01 ·RUL 1 +a 02 ·RUL 2 ,
wherein ,RULi Is the remaining service life RUL of the device in state i; variable t stay,i Is the length of time that has remained in state i so far; d (D) i Is the duration of state i; a, a ij The probability of transition from state i to state j is obtained from the state transition matrix; at state 3, RUL is 0; in case of state 2, there are two possible evolution paths, including with probability a 22 In state 2 sustain D 2 -t stay,2 In units of time or a 23 The probability of (2) transitioning directly to state 3; in the first case, RUL is D 2 -t stay,2 The method comprises the steps of carrying out a first treatment on the surface of the In the second case, RUL is 0; therefore, RUL is a 22 ·(D 2 -t stay,2 ) The method comprises the steps of carrying out a first treatment on the surface of the Corresponding calculations obtain RULs for the device in other states.
Compared with the prior art, the fault prediction method for the air conditioning unit equipment based on the HSMM has the following technical effects:
the invention designs a fault prediction method of air conditioning unit equipment based on HSMM, aiming at the fault diagnosis problem of the air conditioning unit equipment, firstly adopting a neural network, filtering and other methods to extract fault characteristics, and combining a principal component analysis method to obtain a principal fault characteristic matrix; then, establishing initial semi-hidden Markov models corresponding to faults of the air conditioning unit respectively, and training by using a Gibbs sampling method to obtain the semi-hidden Markov models; then, a scaling Viterbi algorithm is applied to obtain state estimation of a semi-hidden Markov model corresponding to each device of the air conditioning unit; then, a backward recursion algorithm is applied to synthesize the semi-hidden Markov models of different faults and the semi-hidden Markov models of the air conditioning unit and equipment thereof, so that the calculated amount is reduced; finally, based on the acquisition of the current state estimation value of the semi-hidden Markov model corresponding to the air conditioning unit and the equipment thereof, further adopting a backward recursion algorithm to calculate and acquire the residual service life of the air conditioning unit and the equipment thereof; therefore, accurate prediction can be realized for the faults of the air conditioning unit equipment, and the actual working efficiency of the air conditioning unit is ensured.
Drawings
FIG. 1 is a flow chart of a fault prediction method of air conditioning unit equipment based on HSMM designed in the invention;
FIG. 2 is a schematic diagram of a semi-hidden Markov model of an air supply fan constructed in accordance with the present invention;
FIG. 3 is a schematic representation of the state transition process established in the design of the present invention corresponding to the combined state of cooling coil state 3;
FIG. 4 is a schematic illustration of the state estimation and actual state of an air conditioning unit synthesized in the design of the present invention;
FIG. 5 is a schematic representation of a predictive air conditioning unit RUL in accordance with the present invention.
Detailed Description
The following describes the embodiments of the present invention in further detail with reference to the drawings.
The invention designs a fault prediction method of an air conditioning unit device based on a dynamic hidden Markov model (Dynamic hidden Markov model, DHMM) and discrete statistical process control (Discrete statistical process control, DSPC), which is implemented in actual application as shown in fig. 1, and specifically comprises the following steps A to D.
Step A, based on the time sequence changes of the related variables of the faults of the air conditioning unit, constructing time sequence changes of the fault characteristics corresponding to the cooling coil and the fan on the air conditioning unit, and further obtaining a fault characteristic matrix X corresponding to the cooling coil and the fan cc 、X sf And further aims at the fault characteristic matrix X cc 、X sf Respectively executing fault characteristic dimension reduction operation to obtain main fault characteristic matrixes O corresponding to the cooling coils and the fans respectively cc 、O sf Step B is then entered.
In practical applications, the step a is specifically performed as follows steps A1 to A5.
Step A1. Collecting the flow of refrigerating water of an air conditioning unit
Figure GDA0004143077470000111
Air flow rate->
Figure GDA0004143077470000112
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA0004143077470000113
And->
Figure GDA0004143077470000114
Respectively used as the time sequence change of the fault characteristics on the air conditioning unit.
At the same time, the air supply temperature T of the air conditioning unit is collected a,sup And a preset value T thereof sup,spt The difference corresponds to the time sequence change delta T from time 1 to time K sup,spt (1)、…、ΔT sup,spt (K) As a fault characteristic time sequence change on the air conditioning unit, then step A2 is entered.
A2, constructing the refrigerating water flow corresponding to the training cooling coil
Figure GDA0004143077470000121
Flow of mixed wind->
Figure GDA0004143077470000122
Temperature T of mixed air a,mix Delivery temperature T with chilled water chw,sup The difference is used as input, and the estimated value of the heat exchange quantity of the cooling coil is +.>
Figure GDA0004143077470000123
Reverse-transfer neural network for output and application +.>
Figure GDA0004143077470000124
T a,mix -T chw,sup Corresponding to time sequence changes from time 1 to time K respectively, obtaining a heat exchange quantity estimated value +.>
Figure GDA0004143077470000125
Corresponding to the time sequence change from time 1 to time K;
then according to the specific heat capacity c of the chilled water pw And the chilled water return temperature T chw,rn The following formula applies:
Figure GDA0004143077470000126
Obtaining the heat exchange quantity Q cc Corresponding to the time sequence change from time 1 to time K, and then
Figure GDA0004143077470000127
Obtaining DeltaQ cc Time sequence change DeltaQ corresponding to time 1 to time K cc (1)、…、ΔQ cc (K) As a result ofAnd (3) carrying out time sequence change for the fault characteristic on the air conditioning unit, and entering a step A3.
Next, particle filtering or the like is adopted to estimate certain parameters irrelevant to the working mode and the working load of the equipment as fault characteristics, namely, the following step A3 is executed.
A3, obtaining an estimated value of the inner diameter of the upper pipe wall of the cooling coil pipe on the air conditioning unit by applying an ash box model and a particle filtering method
Figure GDA0004143077470000128
And rib area estimation +.>
Figure GDA0004143077470000129
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA00041430774700001210
And->
Figure GDA00041430774700001211
Obtaining an efficiency estimate of the fan +.>
Figure GDA00041430774700001212
Corresponding to the time sequence change from time 1 to time K
Figure GDA00041430774700001213
Step A4 is then entered.
In practical applications, the step A3 is specifically performed as follows steps A3-1 to A3-6.
Step A3-1, constructing a cooling coil corresponding to the inner diameter d of the tube wall tube,in Area A of the rib fin The ash bin model of (2) is as follows:
Figure GDA00041430774700001214
wherein ,Ea,mix ＝(1.006T a,mix )+W a,mix (1.84T a,mix +2501),
E a,sup ＝(1.006T a,sup )+W a,sup (1.84T a,sup +2501),
LMTD＝(ΔT sup -ΔT rn )/(lnΔT sup -lnΔT rn )，
τ＝(A fin +A tube,out )/A tube,in ，
α a,e ＝α a [1-A fin (1-η fin )/(A tube,out +A fin )]，
Figure GDA0004143077470000131
wherein ,Wa,mix and Wa,sup Respectively showing the humidity of the mixed air of the air conditioning unit and the humidity of the cooled air supply, E a,mix and Ea,sup Respectively representing the enthalpy value of mixed air of an air conditioning unit and the enthalpy value delta of cooled air supply tube Represents the thickness of the cooling coil pipe, lambda tube Representing the coefficient of thermal conductivity of the cooling coil pipe, R f Indicating the thermal resistance of the cooling coil due to dust, A tube,in Indicating the inner surface area of the cooling coil pipe, A tube,out Represents the outer surface area, deltaT, of the cooling coil pipe sup ＝T a,sup -T chw,sup ,ΔT rn ＝T a,mix -T chw,rn ，T a,sup T represents the air supply temperature chw,rn Represents the return water temperature of chilled water in the cooling coil, alpha a Representing the convective heat transfer coefficient, eta of the rib fin The fin efficiency is represented, the physical coefficient is represented by A,
Figure GDA0004143077470000132
indicating the cooling water flow rate in the cooling coil, ψ indicating the heat flux; step A3-2 is then entered.
Step A3-2. Inner diameter d of pipe wall tube,in Area A of the rib fin As states in the filter frame, a particle filter frame corresponding to the cooling coil is constructed as follows:
cooling coil state equation: x is x cc (t+1)＝x cc (t)+v cc (t)
wherein ,xcc (t+1) represents the state of the cooling coil at the time t+1, x cc (t) represents the state of the cooling coil corresponding to the time t, and the state of the cooling coil comprises the inner diameter d of the tube wall tube,in Area A of the rib fin ，v cc (t) represents the process noise of the cooling coil corresponding to time t, the process noise comprising the process noise of the inner diameter of the tube wall and the process noise of the fin area;
cooling coil measurement equation: z cc (t)＝h(x cc (t))+w cc (t)
wherein ,
Figure GDA0004143077470000133
Figure GDA0004143077470000134
wherein ,wcc (t) represents the measured noise of the cooling coil corresponding to time t, and then step A3-3 is entered.
Step A3-3, based on the cooling coil state equation and the cooling coil measurement equation, using a particle filtering method to represent the posterior probability of a random event by a group of weighted random samples, namely particles, and obtaining the estimated value of the inner diameter of the upper tube wall of the cooling coil on the air conditioning unit from the sensor readings containing noise in the formula
Figure GDA0004143077470000135
And rib area estimation +.>
Figure GDA0004143077470000136
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure GDA0004143077470000141
And->
Figure GDA0004143077470000142
Step A3-4 is then entered.
Step A3-4, constructing a fan to include the efficiency e sf The ash bin model of (2) is as follows:
Figure GDA0004143077470000143
wherein ,Qsf Representing the power of the fan; ρ air Representing the density of air;
Figure GDA0004143077470000144
a preset value representing the air supply quantity; ΔP sf,des Representing the pressure drop of the fan, c 1 、c 2 、c 3 、c 4 、c 5 Coefficients of each term in the polynomial, which can be found in a manual of fans, or obtained based on data training using a nonlinear least squares method; step A3-5 is then entered.
Step A3-5, constructing a particle filter frame corresponding to the fan as follows:
fan state equation: e, e sf (t+1)＝e sf (t)+v sf (t)
wherein ,esf (t+1) represents the efficiency of the fan at time t+1, v sf (t) represents process noise of the fan corresponding to the time t;
the fan measurement equation is as follows:
Figure GDA0004143077470000145
wherein ,
Figure GDA0004143077470000146
w sf (t) represents the measurement noise of the fan at the time t, Q sf (t) represents the power of the fan at the time corresponding to t+1; step A3-6 is then entered.
Step A3-6. Based on the fan State equation and the fan measurement equation, obtaining the fan corresponding to t+1 from the noisy measurement value and the weighted average of the fan efficiency estimates corresponding to the time t by using a particle filtering methodEfficiency estimate of the etch
Figure GDA0004143077470000147
Further obtain the efficiency estimation of the fan>
Figure GDA0004143077470000148
Time sequence variation corresponding to time 1 to time K>
Figure GDA0004143077470000149
A4, obtaining fault characteristic matrixes X corresponding to the cooling coil and the fan respectively cc 、X sf The method comprises the following steps:
Figure GDA00041430774700001410
/>
Figure GDA0004143077470000151
step A5 is then entered.
Step A5. For the failure feature matrix X cc 、X sf Respectively executing PCA principal component analysis method, and aiming at fault feature dimension reduction operation, obtaining a principal fault feature matrix O corresponding to the cooling coil and the fan respectively cc 、O sf 。
For PCA principal component analysis, the dimension reduction process is to extract principal components in the feature matrix to form a new low-dimension matrix, and the specific steps are as follows:
let the feature matrix be
Figure GDA0004143077470000152
Wherein M is the number of feature variables in the feature matrix, and K is the length of the time sequence. First, each row of the feature matrix X is subjected to a decentralization operation, such as the ith row X i ＝[x 1 (i) x 2 (i) … x M (i)]
Figure GDA0004143077470000153
The covariance matrix of the samples in the matrix is then calculated.
C＝XX T ，
Then, the eigenvalues of the covariance matrix C and the corresponding eigenvectors are found. The eigenvalues are arranged according to the size, eigenvectors corresponding to the first k eigenvalues are taken to form a matrix P, and the k eigenvalues are used as main components and need to contain information of all components with the content of more than 90%. For cooling coils, the first three main components form a matrix O cc Can cover 98.4% of information, thus O cc Is used as an observation of a semi-hidden Markov model of the various faults of the cooling coil. Matrix O composed of first two main components of air supply fan sf 96.7% of the information can be covered and thus be used as an observation of the semi-hidden markov model for different faults of the blower fan.
Step B, according to the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf Constructing initial semi-hidden Markov models corresponding to different faults of the cooling coil respectively, constructing initial semi-hidden Markov models corresponding to different faults of the fan respectively, and based on a main fault characteristic matrix O cc and Osf Training the preset proportion training data of the cooling coil by using a Gibbs sampling method aiming at each initial semi-hidden Markov model respectively to obtain the semi-hidden Markov models corresponding to different faults of the cooling coil and the semi-hidden Markov models corresponding to different faults of the fan respectively, and then entering the step C.
In practical application, the step B includes the following steps B1 to B3.
Step B1 based on the main fault characteristic matrix O corresponding to the cooling coil cc As an observation, four faults corresponding to the cooling coil include the occurrence of a deposit f in the cooling coil pipe tube,cc Dust backlog f on cooling coil rib fin,cc The chilled water valve of the cooling coil is blocked in the full-open statef vlvso,cc The chilled water valve of the cooling coil is blocked in the full-closed state f vlvsc,cc And determining initial semi-hidden Markov models corresponding to different faults of the cooling coil based on the normal, initial fault, middle fault and late fault corresponding to the dynamic fault change, wherein each initial semi-hidden Markov model has four states corresponding to the normal, initial fault, middle fault and late fault respectively.
At the same time, based on the main fault characteristic matrix O corresponding to the fan sf As an observation, three faults corresponding to the fan include a fan efficiency drop f ef,sf The fan being locked at a certain speed f fs,sf The fan is thoroughly damaged f cf,sf Determining initial semi-hidden Markov models corresponding to different faults of the fan based on normal, initial fault, middle fault and late fault corresponding to the dynamic fault change, wherein each initial semi-hidden Markov model has four states corresponding to normal, initial fault, middle fault and late fault respectively; and then proceeds to step B2.
Here, the designs correspond to the number 0,1,2 and 3 respectively for the normal, early, middle and late failure states, so that the semi-hidden markov model of the fan has four states, such as shown in fig. 2, wherein the number 0 represents the normal state, the numbers 1,2 and 3 correspond to three failure grades respectively, and the number 3 is the failure degree indicating that the equipment is at the late failure stage, and the equipment can be considered damaged at this time. Thus, from the current time to the time at which the fault enters the fault level corresponding to number 3, it can be considered the remaining service life of the device at the fault.
B2, respectively constructing a state transition matrix P= { a for each initial semi-hidden Markov model corresponding to different faults of the cooling coil and different faults of the fan based on each semi-hidden Markov model having four states ij },i＝1、…、4,j＝1、…、4，{a ij -representing the probability of the device transitioning from the ith state to the jth state; likelihood function b of observations for different states can be constructed based on observations O i (O (s)), i=1,..4, which indicates the probability that the observed value O(s) belongs to the state i, and then proceeds to step B3.
B3, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf Training the initial semi-hidden Markov models by using a Gibbs sampling method according to the preset proportion training data, such as 2/3 training data, combining the initial semi-hidden Markov models, the state transition matrix and the likelihood function of the observed value, which are respectively corresponding to different faults, so as to obtain the semi-hidden Markov models respectively corresponding to the different faults of the cooling coil and the different faults of the fan, and then entering the step C.
Wherein, it is first assumed that the prior distribution of each parameter in the semi-hidden Markov model is a conjugate distribution, and then the posterior estimation of the parameter is deduced based on Bayesian rules. Then updating the estimated value of the parameter through iteration, in the iteration process, firstly selecting one parameter, under the condition of fixing other parameters, generating a random number as the estimated value of the parameter based on posterior distribution of the parameter, and then sequentially and iteratively updating each parameter until meeting convergence conditions, finally obtaining the estimated value of the parameter of the semi-hidden Markov model, wherein the estimated value comprises 1) initial probability vectors pi of the semi-hidden Markov model corresponding to different faults of equipment respectively i I=1,..where N, N is the number of states, determining the probability of belonging to different states at the initial moment; 2) State transition matrix p= { a ij I=1, …, N; j=1, …, N, representing the probability of transitioning between different states; 3) Likelihood function b of observed value i (O (t)), t=1, …, K, where K is the length of the state sequence and O (t) is the observed value at time t; 4) The distribution parameter of the state duration determines the length of the duration.
In a specific application, in order to estimate parameters of the semi-hidden markov model, a gibbs sampling algorithm is used, and the steps of the algorithm are as follows:
1) The a priori distribution of each parameter is assumed to be a conjugate distribution to ensure that the posterior distribution of the parameter obtained by the iterative updating of gibbs sampling is still of the same type. Such as assuming an initial probabilityThe vector pi obeys the dirichlet distribution, pi-Dir (alpha) 1 ,...,α N), wherein ai Is a distribution parameter;
2) Given training data, a posterior distribution of parameters is obtained according to a bayesian rule, for example, the posterior distribution of an initial probability vector pi is pi| to Dir (n 1 +α 1 ,...,n N +α N ) Wherein parameter n i Is the number of points belonging to state i at the initial time.
3) In the iteration process, a parameter is selected firstly, under the condition of fixing other parameters, a random number is generated based on posterior distribution of the parameter to serve as an estimated value of the parameter, and then each parameter is sequentially and iteratively updated until convergence conditions are met, and finally the estimated value of the parameter of the semi-hidden Markov model is obtained.
Step C, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf And (3) obtaining a semi-hidden Markov model corresponding to different faults of the cooling coil and different faults of the fan respectively and corresponding state estimation from time 1 to time K by adopting a scaling Viterbi algorithm, and then entering the step D.
In practical applications, the step C is to obtain the state estimation of the semi-hidden markov models corresponding to the time 1 to the time K according to the following steps C1 to C8, i.e. obtain the state estimation of the semi-hidden markov models corresponding to the different faults of the cooling coil and the different faults of the fan respectively corresponding to the time 1 to the time K.
Step c1. When the initialization time t=1, the hidden state is the probability maximum value in all possible state transition paths of i e {1, 2, …, N }.
Step C2. forward procedure: the following formula is adopted:
Figure GDA0004143077470000181
calculating the probability that the observation value corresponds to the state i of the equipment from the time t-d+1 to the time tRate delta i (t) wherein,
Figure GDA0004143077470000182
p i (d) Representing the probability that the device is in state i for d time units; b i (O (s)) represents the probability that the observed value O(s) belongs to the state i.
Step C3. is according to the following formula:
Figure GDA0004143077470000183
obtaining
Figure GDA0004143077470000184
wherein />
Figure GDA0004143077470000185
And is based on
Figure GDA0004143077470000186
Calculating a scaling factor->
Figure GDA0004143077470000187
wherein ,bi (O (t)) represents an observation likelihood function corresponding to a time t in the semi-hidden Markov model; when k (t) approaches ≡, the scaling factor k (t) is set to 1 +.>
Figure GDA0004143077470000188
If k (t) does not approach ≡, the value corresponding to k (t) is retained.
And step C4, adjusting various parameters in the semi-hidden Markov model based on the scaling factor k (t).
Step C5. obtains the maximized probability
Figure GDA0004143077470000189
Duration d of (2), where d=1, …, M i ，M i Is the upper limit for the duration of state i.
Step C6. obtains maximized delta i (t)a ij Corresponding toState i of (i), i.e
Figure GDA00041430774700001810
wherein ,aij Representing the probability that the device transitions from the i-th state to the j-th state.
Step C7. backward procedure: calculating the probability of all states at time K and determining the most probable state delta i (K)；
Step C8. obtains the state estimation of the failure semi-hidden Markov model from time 1 to time K by backtracking the recorded d and i
Figure GDA0004143077470000191
And D, integrating the semi-hidden Markov models corresponding to the faults of the equipment respectively by applying a backward recursion algorithm to the cooling coil and the fan to obtain the air conditioning unit, the semi-hidden Markov models corresponding to the cooling coil and the fan respectively, further obtaining the current state estimated values of the air conditioning unit, the cooling coil and the fan corresponding to the semi-hidden Markov models, and then step E.
In practical applications, in the step D, the following steps D1 to D6 are executed for the cooling coil, the fan and the air conditioning unit, respectively, to obtain the air conditioning unit, the semi-hidden markov model corresponding to the cooling coil and the fan, and the current state estimation value.
D1, generating N according to state transition matrix of semi-hidden Markov model corresponding to different faults of equipment by using Monte Carlo simulation smp A group state sequence, according to the state sequences, randomly selecting a value as the value of the fault parameter in the interval of the fault parameter corresponding to the discrete state, substituting the value of the fault parameter into the ash box model of the cooling coil and the ash box model of the fan, calculating the corresponding equipment capacity, wherein the equipment capacity is the state of the equipment, and further obtaining the fault state s f And device state s dev Sampling points of the relationship between the two; associating each device with a fault state s f Together fault co-operating state s js To show that, by means of the method,the characterization of the fault joint state s can be further obtained js And device state s dev Sampling points of the relationship between the two.
Step D2. is to extract two mapping matrices for each device according to the semi-hidden Markov model corresponding to the different faults of the device and the state estimation corresponding to the time 1 to the time K, respectively, from the sampling points obtained in step D1, to embody the fault joint state s js And device state s dev Relationship between them.
Consider, for example, two failures of a cooling coil, including pipe deposit f tube,cc And dust backlog f on the rib fin,cc And a failure of the fan, i.e. a drop in fan efficiency f eff,sf . The three faults respectively correspond to three fault related parameters, namely the inner diameter d of the pipe wall tube,in Surface area of rib A fin And fan efficiency e sf . The joint probability of failure is defined as s js ＝(s tube ,s fin ,s sf), wherein stube 、s fin and ssf Respectively f tube,cc 、f fin,cc F eff,sf Is a fault condition of (a). Thus, the fault joint state depends on different fault states, one for each bit. Wherein, "wherein indicates that the fault did not occur," indicates that the, "" the indicates the extent of the fault. For example, s js =1 corresponds to "01", meaning f eff,sf Is failure level 1, and f tube,cc 、f fin,cc No occurrence occurs.
The fault joint state s is embodied js And device state s dev Two mapping matrices of the relationship between each other, the first mapping matrix being represented by a given s dev Time s js Probability of each state, the matrix is defined as M js|dev ＝{m ij,js|dev}, wherein ,
m ij,js|dev ＝Pr(s js ＝j|s dev ＝i)＝n ij,dev /n i,dev 。
in the formula,sdev Is the status of the device; s is(s) js Is the joint state of the fault corresponding to the device；n ij,dev Is corresponding to s dev =i and s js Number of sampling points of =j, n i,dev Is corresponding to state s dev Number of points=i.
Corresponding second mapping matrix M dev|js ＝{m ji,dev|js The expression of a given s js Time s dev Probability of each state, where m ji,dev|js ＝Pr(s dev ＝i|s js ＝j)。
And step D3, integrating parameters of the semi-hidden Markov models corresponding to the faults respectively according to the two mapping matrixes of each device to obtain the semi-hidden Markov model corresponding to the device.
Step D4. applies a backward recursion algorithm according to the distribution obeyed by the state durations corresponding to the different faults respectively, and obtains the duration of each state of each device under the faults respectively by calculating the evolution time among different joint states.
For example, the progression of the degree of failure over time is always changing in a more severe direction. The digital sum of the fault joint states also evolves from a small to a large direction. For example, as shown in FIG. 3, when the cooling coil is in state 3, its two failure states s tube,cc ,s fin,cc May be { s } tube,cc ,s fin,cc }＝{2,1},{s tube,cc ,s fin,cc }＝{2,2},{s tube,cc ,s fin,cc }＝{3,1},or{s tube,cc ,s fin,cc = {3,2}. I.e. s cc =3 corresponds to the joint state s js,cc = { '21, "22,"31, "32' }. Considering that the digital sum of the joint states evolves from a small to large direction, then from the joint state s js,cc = '21' evolves to a joint state s js,cc The time of= '32' can be considered as cooling coil state s cc Expected duration of=3. To calculate the slave joint state
Figure GDA0004143077470000201
To the joint state->
Figure GDA0004143077470000202
Is required to know the evolution time from the joint state +.>
Figure GDA0004143077470000203
To the joint state->
Figure GDA0004143077470000204
Is a time of evolution of (2). Thus, the backward recursion algorithm is used to obtain the duration of the device state by calculating the evolution time between the different joint states, and the execution of the specific step D4 is as follows:
in step D4, the joint state j belongs to the state k of the device, and when the digital sum of the joint state j is D max,k ：
Figure GDA0004143077470000205
wherein dmax,k Is the maximum value of the sum of the digits of the joint state corresponding to the equipment state k;
Figure GDA0004143077470000206
is the set of joint states where the sum of the digits of the joint states is the maximum.
When the digital sum of the joint state j is d max,k -1：
Figure GDA0004143077470000207
…
When the digital sum of the joint state j is d min,k ：
Figure GDA0004143077470000211
wherein Dj,js Is the duration of the joint state j;
Figure GDA0004143077470000212
is a slave unitThe combination state j to the combination state->
Figure GDA0004143077470000213
Is a time of evolution of (2); when->
Figure GDA0004143077470000214
When (I)>
Figure GDA0004143077470000215
The duration of time that is considered to be the device state k; will->
Figure GDA0004143077470000216
The number of joint states of (a) is defined as n min,k The method comprises the steps of carrying out a first treatment on the surface of the If n min,k =1, the duration of device state k is +.>
Figure GDA0004143077470000217
If n min,k >1, then the average value of the durations corresponding to the different joint states j>
Figure GDA0004143077470000218
May be considered as the duration of state k.
Step D5. according to embodying the Fault Joint State s js And device state s dev The two mapping matrixes of the relation between the two mapping matrixes, the duration time of each state of the equipment under each fault and the transition matrix of the joint state respectively generate a state transition matrix of a semi-hidden Markov model corresponding to the equipment, and the state transition matrix is obtained according to the following formula for the equipment;
Figure GDA0004143077470000219
wherein ,sjs Is the joint state of the devices; p(s) dev (t+1)＝j|s js (t+1) =n) and P(s) js (t)＝m|s dev (t) =i) from the mapping matrix.
Step D6. obtains the current state estimation value of the semi-hidden Markov model corresponding to the device according to the states of different faults of the semi-hidden Markov model corresponding to the device.
Wherein, based on the three gradual faults of the device being independent, pr(s) js (t)＝i):
Figure GDA00041430774700002110
wherein ,fl ,l＝1,...,N dev Is N of the device dev A fault; joint state s js (t) =i corresponds to N dev Failure f l ,l＝1,...,N dev In state i l ,l＝1,...,N dev In obtaining the joint state s js After the probability of each state, the state of the device at time t is estimated as follows
Figure GDA00041430774700002111
Based on the mapping matrix m ij,dev|js To estimate the number of times the sample is to be processed,
Figure GDA00041430774700002112
wherein Njs Is the number of joint states, and as shown in fig. 4, the current state estimation value of the device is obtained by synthesizing state estimation values of different faults.
And E, calculating to obtain the residual service lives of the air conditioning unit, the cooling coil and the fan by adopting a backward recursion algorithm according to the air conditioning unit, the semi-hidden Markov model corresponding to the cooling coil and the fan and the current state estimated value of each device.
In practical application, the specific steps of the step E are as follows: based on the parameters of the dynamic hidden Markov model and the estimated values of the states, a backward recursion algorithm is employed to estimate the remaining useful life. The algorithm estimates the parameters and states of the semi-hidden Markov model corresponding to each device of the air conditioning unit based on the fifth step, and adopts a recursion method to estimate the time length from the current moment to the moment when the device is in fault, namely the residual service life of each device of the air conditioning unit.
The HSMM of each device has four states, with state 0 corresponding to the normal state and states 1, 2, 3 corresponding to three levels of failure. These fault levels are determined by the severity of all of their faults. For example, the upper limit of the degree of failure of the blower fan is level 3 (state 3), that is, when the failure state is 3, the device is considered to be not operating normally. Thus, the remaining useful life of the device in state 3 is 0; when the fault is in state 2, it is possible for the fault to stay in its current state 2 or transition to state 3; similarly, when a fault is in state 1, the fault may stay in its current state 1, or other states, including state 2 and state 3. The method is based on four states of the semi-hidden Markov model corresponding to the fans respectively, and comprises the following specific steps:
state 3 RUL 3 ＝0,
State 2 RUL 2 ＝a 22 ·(D 2 -t stay,2 ),
State 1 RUL 1 ＝a 11 ·(D 1 -t stay,1 +RUL 2 )+a 12 ·RUL 2 ,
State 0 RUL 0 ＝a 00 ·(D 0 -t stay,0 +RUL 1 )+a 01 ·RUL 1 +a 02 ·RUL 2 ,
wherein ,RULi Is the remaining service life RUL of the device in state i; variable t stay,i Is the length of time that has remained in state i so far; d (D) i Is the duration of state i; a, a ij The probability of transition from state i to state j is obtained from the state transition matrix; at state 3, RUL is 0; in case of state 2, there are two possible evolution paths, including with probability a 22 In state 2 sustain D 2 -t stay,2 In units of time or a 23 The probability of (2) transitioning directly to state 3; in the first case, RUL is D 2 -t stay,2 The method comprises the steps of carrying out a first treatment on the surface of the In the second case, RUL is 0; therefore, RUL is a 22 ·(D 2 -t stay,2 ) The method comprises the steps of carrying out a first treatment on the surface of the Corresponding calculations obtain RULs for the device in other states.
The fault prediction method of the air conditioning unit equipment based on the HSMM is applied to practice, such as designing and collecting 45793 samples in the No. 11-22 of a small building from 6 months, wherein half of data are used as training data, and the other half of data are used as test data. Parameters of a semi-hidden Markov model of the cooling coil were trained using the Gibbs sampling method. To measure the error magnitude, F-measure is used. This is because F-measure depends not only on the rate of missing report (failure state misdiagnosis is normal state) but also on the rate of false report (one failure state is misreported as normal state or other failure state), and the closer to 1, the smaller the rate of missing report and the rate of false report. The scaling mechanism in the patent realizes dynamic adjustment of the scaling factor through a mechanism for switching the scaling factor under different conditions, avoids the condition that the forward probability after scaling approaches to +.. The state estimates for F-measure are 0.998,0.946,0.835 and 0.928, respectively, which are close to 1 and therefore relatively accurate. By adopting the synthesis method provided by the invention, the state estimation of different faults can be synthesized to directly obtain the state estimation of the air conditioning unit. As shown in fig. 4, the error between the state estimation of the combined air conditioning unit and its actual state is small. The remaining useful life may be estimated based on the estimated value of the current state and parameters of the semi-hidden Markov model. As shown in fig. 5, the remaining useful life of the cooling coil is estimated. Wherein the relative accuracy of the predictions is [0.662,0.595,0.845], which correspond to the front data, middle data and end data, respectively. These 3 values are close to 1 and therefore the prediction accuracy is relatively high.
The embodiments of the present invention have been described in detail with reference to the drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the spirit of the present invention.

Claims (8)

1. The fault prediction method of the air conditioning unit equipment based on the HSMM is characterized by comprising the following steps of:
step A, based on the time sequence changes of the related variables of the faults of the air conditioning unit, constructing time sequence changes of the fault characteristics corresponding to the cooling coil and the fan on the air conditioning unit, and further obtaining a fault characteristic matrix X corresponding to the cooling coil and the fan cc 、X sf And further aims at the fault characteristic matrix X cc 、X sf Respectively executing fault characteristic dimension reduction operation to obtain main fault characteristic matrixes O corresponding to the cooling coils and the fans respectively cc 、O sf Then enter step B;
step B, according to the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf Constructing initial semi-hidden Markov models corresponding to different faults of the cooling coil respectively, constructing initial semi-hidden Markov models corresponding to different faults of the fan respectively, and based on a main fault characteristic matrix O cc and Osf Training the preset proportion training data of the cooling coil by using a Gibbs sampling method aiming at each initial semi-hidden Markov model respectively to obtain semi-hidden Markov models corresponding to different faults of the cooling coil and semi-hidden Markov models corresponding to different faults of the fan respectively, and then entering the step C;
step C, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf The method comprises the steps of (1) obtaining a semi-hidden Markov model corresponding to different faults of a cooling coil and different faults of a fan respectively and corresponding to state estimation from time 1 to time K by adopting a scaling Viterbi algorithm, and then entering a step D;
step D, integrating the semi-hidden Markov models corresponding to the faults of the equipment respectively by applying a backward recursion algorithm to the cooling coil and the fan to obtain an air conditioning unit, the semi-hidden Markov models corresponding to the cooling coil and the fan respectively, further obtaining the current state estimation values of the air conditioning unit, the cooling coil and the fan corresponding to the semi-hidden Markov models, and then entering the step E;
and E, calculating to obtain the residual service lives of the air conditioning unit, the cooling coil and the fan by adopting a backward recursion algorithm according to the air conditioning unit, the semi-hidden Markov model corresponding to the cooling coil and the fan and the current state estimated value of each device.
2. The method for predicting faults of air conditioning unit equipment based on the HSMM according to claim 1, wherein the step a includes the following steps A1 to A5:
step A1. Collecting the flow of refrigerating water of an air conditioning unit
Figure FDA0004143077460000011
Air flow rate->
Figure FDA0004143077460000012
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure FDA0004143077460000013
And->
Figure FDA0004143077460000014
Respectively used as fault characteristic time sequence changes on the air conditioning unit;
at the same time, the air supply temperature T of the air conditioning unit is collected a,sup And a preset value T thereof sup,spt The difference corresponds to the time sequence change delta T from time 1 to time K sup,spt (1)、…、ΔT sup,spt (K) The method comprises the steps of taking the fault characteristic time sequence change as a fault characteristic time sequence change on an air conditioning unit, and then entering a step A2;
a2, constructing the refrigerating water flow corresponding to the training cooling coil
Figure FDA0004143077460000021
Flow of mixed wind->
Figure FDA0004143077460000022
Temperature T of mixed air a,mix And freezingTemperature T of water delivery chw,sup The difference is used as input, and the estimated value of the heat exchange quantity of the cooling coil is +.>
Figure FDA0004143077460000023
Reverse-transfer neural network for output and application +.>
Figure FDA0004143077460000024
T a,mix -T chw,sup Corresponding to time sequence changes from time 1 to time K respectively, obtaining a heat exchange quantity estimated value +.>
Figure FDA0004143077460000025
Corresponding to the time sequence change from time 1 to time K;
then according to the specific heat capacity c of the chilled water pw And the chilled water return temperature T chw,rn The following formula applies:
Figure FDA0004143077460000026
obtaining the heat exchange quantity Q cc Corresponding to the time sequence change from time 1 to time K, and then
Figure FDA0004143077460000027
Obtaining DeltaQ cc Time sequence change DeltaQ corresponding to time 1 to time K cc (1)、…、ΔQ cc (K) As the time sequence change of the fault characteristic on the air conditioning unit, entering into the step A3;
a3, obtaining an estimated value of the inner diameter of the upper pipe wall of the cooling coil pipe on the air conditioning unit by applying an ash box model and a particle filtering method
Figure FDA0004143077460000028
And rib area estimation +.>
Figure FDA0004143077460000029
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure FDA00041430774600000210
And->
Figure FDA00041430774600000211
Obtaining an efficiency estimate of the fan +.>
Figure FDA00041430774600000212
Corresponding to the time sequence change from time 1 to time K
Figure FDA00041430774600000213
Then enter step A4;
a4, obtaining fault characteristic matrixes X corresponding to the cooling coil and the fan respectively cc 、X sf The method comprises the following steps:
Figure FDA00041430774600000214
Figure FDA00041430774600000215
then enter step A5;
step A5. For the failure feature matrix X cc 、X sf Respectively executing PCA principal component analysis method, and aiming at fault feature dimension reduction operation, obtaining a principal fault feature matrix O corresponding to the cooling coil and the fan respectively cc 、O sf 。
3. The method for predicting faults of air conditioning unit equipment based on the HSMM according to claim 2, wherein the step A3 comprises the following steps A3-1 to A3-6:
step A3-1, constructing a cooling coil corresponding to the inner diameter d of the tube wall tube,in Area A of the rib fin The ash bin model of (2) is as follows:
Figure FDA0004143077460000031
wherein ,Ea,mix ＝(1.006T a,mix )+W a,mix (1.84T a,mix +2501),
E a,sup ＝(1.006T a,sup )+W a,sup (1.84T a,sup +2501),
LMTD＝(ΔT sup -ΔT rn )/(lnΔT sup -lnΔT rn )，
τ＝(A fin +A tube,out )/A tube,in ，
α a,e ＝α a [1-A fin (1-η fin )/(A tube,out +A fin )]，
Figure FDA0004143077460000032
wherein ,Wa,mix and Wa,sup Respectively showing the humidity of the mixed air of the air conditioning unit, the humidity of the cooled air supply and E a,mix and Ea,sup Respectively representing the enthalpy value of mixed air of an air conditioning unit and the enthalpy value delta of cooled air supply tube Represents the thickness of the cooling coil pipe, lambda tube Representing the coefficient of thermal conductivity of the cooling coil pipe, R f Indicating the thermal resistance of the cooling coil due to dust, A tube,in Indicating the inner surface area of the cooling coil pipe, A tube,out Represents the outer surface area, deltaT, of the cooling coil pipe sup ＝T a,sup -T chw,sup ,ΔT rn ＝T a,mix -T chw,rn ，T a,sup T represents the air supply temperature chw,rn Represents the return water temperature of chilled water in the cooling coil, alpha a Representing the convective heat transfer coefficient, eta of the rib fin The fin efficiency is represented, the physical coefficient is represented by A,
Figure FDA0004143077460000033
indicating the cold in the cooling coilThe flow rate of the frozen water, psi, represents the heat flow density; then enter step A3-2;
step A3-2. Inner diameter d of pipe wall tube,in Area A of the rib fin As states in the filter frame, a particle filter frame corresponding to the cooling coil is constructed as follows:
cooling coil state equation: x is x cc (t+1)＝x cc (t)+v cc (t)
wherein ,xcc (t+1) represents a state of the cooling coil at time t+1, and x represents a state of the cooling coil at time t+1 cc (t) represents the state of the cooling coil corresponding to the time t, and the state of the cooling coil comprises the inner diameter d of the tube wall tube,in Area A of the rib fin ，v cc (t) represents the process noise of the cooling coil corresponding to time t, the process noise comprising the process noise of the inner diameter of the tube wall and the process noise of the fin area;
Cooling coil measurement equation: z cc (t)＝h(x cc (t))+w cc (t)
wherein ,
Figure FDA0004143077460000041
Figure FDA0004143077460000042
wherein ,wcc (t) representing the measured noise of the cooling coil corresponding to the time t, and then entering the step A3-3;
step A3-3, based on the cooling coil state equation and the cooling coil measurement equation, using a particle filtering method to represent the posterior probability of a random event by a group of weighted random samples, namely particles, and obtaining the estimated value of the inner diameter of the upper tube wall of the cooling coil on the air conditioning unit from the sensor readings containing noise in the formula
Figure FDA0004143077460000043
And rib area estimation +.>
Figure FDA0004143077460000044
Time sequence variation +_corresponding to time 1 to time K respectively>
Figure FDA0004143077460000045
And->
Figure FDA0004143077460000046
Then enter step A3-4;
step A3-4, constructing a fan to include the efficiency e sf The ash bin model of (2) is as follows:
Figure FDA0004143077460000047
wherein ,Qsf Representing the power of the fan; ρ air Representing the density of air;
Figure FDA0004143077460000048
a preset value representing the air supply quantity; ΔP sf,des Representing the pressure drop of the fan, c 1 、c 2 、c 3 、c 4 、c 5 Coefficients of each term in the polynomial, which coefficients are found in a manual of fans, or are obtained based on data training using a nonlinear least squares method; then enter step A3-5;
step A3-5, constructing a particle filter frame corresponding to the fan as follows:
fan state equation: e, e sf (t+1)＝e sf (t)+v sf (t)
wherein ,esf (t+1) represents the efficiency of the fan at time t+1, v sf (t) represents process noise of the fan corresponding to the time t;
the fan measurement equation is as follows:
Figure FDA0004143077460000049
wherein ,
Figure FDA00041430774600000410
w sf (t) represents the measurement noise of the fan at the time t, Q sf (t) represents the power of the fan at the time corresponding to t+1; then enter step A3-6;
step A3-6, based on the fan state equation and the fan measurement equation, obtaining an efficiency estimation value of the fan corresponding to the time t+1 from the noise measurement value and a weighted average of the fan efficiency estimation corresponding to the time t by adopting a particle filtering method
Figure FDA0004143077460000051
Further obtain the efficiency estimation of the fan>
Figure FDA0004143077460000052
Time sequence variation corresponding to time 1 to time K>
Figure FDA0004143077460000053
/>
4. The fault prediction method for an air conditioning unit device based on an HSMM according to claim 1, wherein the step B includes the following steps B1 to B3;
step B1 based on the main fault characteristic matrix O corresponding to the cooling coil cc As observation, determining initial semi-hidden Markov models corresponding to different faults of the cooling coil, wherein each initial semi-hidden Markov model has four states and corresponds to normal, initial fault, middle fault and late fault respectively;
at the same time, based on the main fault characteristic matrix O corresponding to the fan sf As observation, determining initial semi-hidden Markov models corresponding to different faults of the fan respectively, wherein each initial semi-hidden Markov model has four states and corresponds to normal, initial fault, middle fault and late fault respectively; then enter step B2;
Step B2, aiming at different faults and fans of the cooling coilInitial semi-hidden Markov models corresponding to different faults respectively, and constructing a state transition matrix P= { a based on four states of each semi-hidden Markov model aiming at each initial semi-hidden Markov model ij },i＝1、…、4,j＝1、...、4，{a ij -representing the probability of the device transitioning from the ith state to the jth state; constructing likelihood function b for observations of different states based on observations O i (O (s)), i=1,..4, indicating the probability that the observed value O(s) belongs to the state i, then proceeding to step B3;
b3, based on the main fault characteristic matrix O corresponding to the cooling coil and the fan respectively cc and Osf And (3) training the initial semi-hidden Markov models by using a Gibbs sampling method according to the preset proportion training data of the cooling coil, combining the initial semi-hidden Markov models, the state transition matrix and the likelihood functions of the observed values, which correspond to different faults respectively, so as to obtain the semi-hidden Markov models corresponding to the different faults of the cooling coil and the different faults of the fan respectively, and then entering the step C.
5. The method for predicting failure of an HSMM based air conditioning unit according to claim 1 or 4, wherein the cooling coil differential failure includes occurrence of a deposit f in the cooling coil tubing tube,cc Dust backlog f on cooling coil rib fin,cc The chilled water valve of the cooling coil is blocked in the full-open state f vlvso,cc The chilled water valve of the cooling coil is blocked in the full-closed state f vlvsc,cc ；
The different faults of the fan comprise the reduction f of the efficiency of the fan ef,sf The fan being locked at a certain speed f fs,sf The fan is thoroughly damaged f cf,sf 。
6. The fault prediction method for air conditioning unit equipment based on the HSMM according to claim 1, wherein the fault prediction method comprises the following steps: in the step C, for each semi-hidden markov model corresponding to different faults of the cooling coil and different faults of the fan, according to the following steps C1 to C8, state estimation from time 1 to time K corresponding to the semi-hidden markov model is obtained, that is, state estimation from time 1 to time K corresponding to the semi-hidden markov model corresponding to different faults of the cooling coil and different faults of the fan is obtained;
step c1, when the initialization time t=1, the hidden state is the probability maximum value in all possible state transition paths of i e {1, 2, …, N };
step C2. forward procedure: the following formula is adopted:
Figure FDA0004143077460000061
calculating the probability delta that the observed value corresponds to the state i of the equipment from the time t-d+1 to the time t i (t) wherein,
Figure FDA0004143077460000062
p i (d) Representing the probability that the device is in state i for d time units; b i (O (s)) represents the probability that the observed value O(s) belongs to the state i;
step C3. is according to the following formula:
Figure FDA0004143077460000063
/>
obtaining
Figure FDA0004143077460000064
wherein />
Figure FDA0004143077460000065
And based on->
Figure FDA0004143077460000066
Calculating a scaling factor->
Figure FDA0004143077460000067
wherein ,bi (O (t)) represents that the observed value of the corresponding time t in the semi-hidden Markov model is similarA natural function; when k (t) approaches ≡, the scaling factor k (t) is set to 1 +.>
Figure FDA0004143077460000068
If k (t) is not approaching ≡, reserving the value corresponding to k (t);
step C4, adjusting various parameters in the semi-hidden Markov model based on a scaling factor k (t);
step C5. obtains the maximized probability
Figure FDA0004143077460000069
Duration d of (2), where d=1, …, M i ，M i Is the upper limit of the duration of state i;
step C6. obtains maximized delta i (t)a ij Corresponding state i, i.e
Figure FDA0004143077460000071
wherein ,aij Representing a probability that the device transitions from the ith state to the jth state;
step C7. backward procedure: calculating the probability of all states at time K and determining the most probable state delta i (K)；
Step C8. obtains the state estimation of the failure semi-hidden Markov model from time 1 to time K by backtracking the recorded d and i
Figure FDA0004143077460000072
7. The fault prediction method for air conditioning unit equipment based on the HSMM according to claim 1, wherein the fault prediction method comprises the following steps: in the step D, the following steps D1 to D6 are executed for the cooling coil, the fan and the air conditioning unit, respectively, to obtain the air conditioning unit, the semi-hidden markov model corresponding to the cooling coil and the fan thereof, and the current state estimation value;
D1. using Monte Carlo simulation to correspond to semi-hidden Markov models based on different faults of the deviceState transition matrix, generating N smp A group state sequence, according to the state sequences, randomly selecting a value as the value of the fault parameter in the interval of the fault parameter corresponding to the discrete state, substituting the value of the fault parameter into the ash box model of the cooling coil and the ash box model of the fan, calculating the corresponding equipment capacity, wherein the equipment capacity is the state of the equipment, and further obtaining the fault state s f And device state s dev Sampling points of the relationship between the two; associating each device with a fault state s f Together fault co-operating state s js To indicate that the characterization of the fault joint state s can be further obtained js And device state s dev Sampling points of the relationship between the two;
step D2. is to extract two mapping matrices for each device according to the semi-hidden Markov model corresponding to the different faults of the device and the state estimation corresponding to the time 1 to the time K, respectively, from the sampling points obtained in step D1, to embody the fault joint state s js And device state s dev A relationship between;
the fault joint state s is embodied js And device state s dev Two mapping matrices of the relationship between each other, the first mapping matrix being represented by a given s dev Time s js Probability of each state, the matrix is defined as M js|dev ＝{m ij,js|dev}, wherein ,
m ij,js|dev ＝Pr(s js ＝j|s dev ＝i)＝n ij,dev /n i,dev ；
in the formula,sdev Is the status of the device; s is(s) js Is the joint state of the fault to which the device corresponds; n is n ij,dev Is corresponding to s dev =i and s js Number of sampling points of =j, n i,dev Is corresponding to state s dev Number of points=i;
corresponding second mapping matrix M dev|js ＝{m ji,dev|js The expression of a given s js Time s dev Probability of each state, where m ji,dev|js ＝Pr(s dev ＝i|s js ＝j)；
Step D3, integrating parameters of the semi-hidden Markov models corresponding to the faults respectively according to the two mapping matrixes of each device to obtain the semi-hidden Markov model corresponding to the device;
step D4., applying a backward recursion algorithm according to the distribution obeyed by the state durations corresponding to different faults respectively, and obtaining the duration of each state of each device under each fault respectively by calculating the evolution time among different joint states; in step D4, the joint state j belongs to the state k of the device, and when the digital sum of the joint state j is D max,k ：
Figure FDA0004143077460000081
wherein dmax,k Is the maximum value of the sum of the digits of the joint state corresponding to the equipment state k;
Figure FDA0004143077460000082
is a set of joint states in which the sum of digits of the joint states is a maximum;
when the digital sum of the joint state j is d max,k -1：
Figure FDA0004143077460000083
When the digital sum of the joint state j is d min,k ：
Figure FDA0004143077460000084
wherein Dj,js Is the duration of the joint state j;
Figure FDA0004143077460000085
from the joint state j to the joint state +.>
Figure FDA0004143077460000086
Is a time of evolution of (2); when->
Figure FDA0004143077460000087
When (I)>
Figure FDA0004143077460000088
The duration of time that is considered to be the device state k; will->
Figure FDA0004143077460000089
The number of joint states of (a) is defined as n min,k The method comprises the steps of carrying out a first treatment on the surface of the If n min,k =1, the duration of device state k is +.>
Figure FDA00041430774600000810
If n min,k >1, then the average value of the durations corresponding to the different joint states j>
Figure FDA00041430774600000811
The duration of what is considered to be state k;
step D5. according to embodying the Fault Joint State s js And device state s dev The two mapping matrixes of the relation between the two mapping matrixes, the duration time of each state of the equipment under each fault and the transition matrix of the joint state respectively generate a state transition matrix of a semi-hidden Markov model corresponding to the equipment, and the state transition matrix is obtained according to the following formula for the equipment;
Figure FDA00041430774600000812
wherein ,sjs Is the joint state of the devices; p(s) dev (t+1)＝j|s js (t+1) =n) and P(s) js (t)＝m|s dev (t) =i) from the mapping matrix;
step D6., obtaining the current state estimation value of the semi-hidden Markov model corresponding to the equipment according to the states of different faults of the semi-hidden Markov model corresponding to the equipment;
wherein, if three gradual faults based on the equipment are independent, pr (s js (t)＝i):
Figure FDA0004143077460000091
wherein ,fl ,l＝1,...,N dev Is N of the device dev A fault; joint state s js (t) =i corresponds to N dev Failure f l ,l＝1,...,N dev In state i l ,l＝1,...,N dev In obtaining the joint state s js After the probability of each state, the state of the device at time t is estimated as follows
Figure FDA0004143077460000092
Based on the mapping matrix m ij,dev|js To estimate the number of times the sample is to be processed,
Figure FDA0004143077460000093
wherein Njs Is the number of joint states, and then obtains the current state estimation value of the equipment based on the states of different faults.
8. The fault prediction method for air conditioning unit equipment based on the HSMM according to claim 1, wherein the fault prediction method comprises the following steps: in the step E, according to the air conditioning unit, the semi-hidden Markov model corresponding to the cooling coil and the fan of the air conditioning unit and the corresponding current state estimation value, a backward recursion algorithm is adopted, and the remaining service lives of the air conditioning unit, the cooling coil and the fan of the air conditioning unit are obtained through calculation in the following manner;
the air conditioning unit and the semi-hidden Markov model corresponding to the air conditioning unit and the equipment respectively have four states, and the method is implemented as follows, wherein state 0 corresponds to a normal state, and states 1, 2 and 3 correspond to three fault levels respectively;
state 3 RUL 3 ＝0,
State 2 RUL 2 ＝a 22 ·(D 2 -t stay,2 ),
State 1 RUL 1 ＝a 11 ·(D 1 -t stay,1 +RUL 2 )+a 12 ·RUL 2 ,
State 0 RUL 0 ＝a 00 ·(D 0 -t stay,0 +RUL 1 )+a 01 ·RUL 1 +a 02 ·RUL 2 ,
wherein ,RULi Is the remaining service life RUL of the device in state i; variable t stay,i Is the length of time that has remained in state i so far; d (D) i Is the duration of state i; a, a ij The probability of transition from state i to state j is obtained from the state transition matrix; at state 3, RUL is 0; in case of state 2, there are two possible evolution paths, including with probability a 22 In state 2 sustain D 2 -t stay,2 In units of time or a 23 The probability of (2) transitioning directly to state 3; in the first case, RUL is D 2 -t stay,2 The method comprises the steps of carrying out a first treatment on the surface of the In the second case, RUL is 0; therefore, RUL is a 22 ·(D 2 -t stay,2 ) The method comprises the steps of carrying out a first treatment on the surface of the Corresponding calculations obtain RULs for the device in other states.
CN202110262737.0A 2021-03-11 2021-03-11 Fault prediction method of air conditioning unit equipment based on HSMM Active CN113158546B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110262737.0A CN113158546B (en) 2021-03-11 2021-03-11 Fault prediction method of air conditioning unit equipment based on HSMM

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110262737.0A CN113158546B (en) 2021-03-11 2021-03-11 Fault prediction method of air conditioning unit equipment based on HSMM

Publications (2)

Publication Number Publication Date
CN113158546A CN113158546A (en) 2021-07-23
CN113158546B true CN113158546B (en) 2023-05-09

Family

ID=76887006

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110262737.0A Active CN113158546B (en) 2021-03-11 2021-03-11 Fault prediction method of air conditioning unit equipment based on HSMM

Country Status (1)

Country Link
CN (1) CN113158546B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104504296A (en) * 2015-01-16 2015-04-08 湖南科技大学 Gaussian mixture hidden Markov model and regression analysis remaining life prediction method
CN109344518A (en) * 2018-10-14 2019-02-15 毛述春 A kind of method for diagnosing faults of base station heat management system
CN109829468A (en) * 2018-04-16 2019-05-31 南京航空航天大学 Civil aircraft Fault Diagnosis of Complex System method based on Bayesian network
CN111981635A (en) * 2020-07-21 2020-11-24 沈阳安新自动化控制有限公司 Central air conditioner fault prediction and diagnosis method adopting double intelligent algorithms

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7904187B2 (en) * 1999-02-01 2011-03-08 Hoffberg Steven M Internet appliance system and method

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104504296A (en) * 2015-01-16 2015-04-08 湖南科技大学 Gaussian mixture hidden Markov model and regression analysis remaining life prediction method
CN109829468A (en) * 2018-04-16 2019-05-31 南京航空航天大学 Civil aircraft Fault Diagnosis of Complex System method based on Bayesian network
CN109344518A (en) * 2018-10-14 2019-02-15 毛述春 A kind of method for diagnosing faults of base station heat management system
CN111981635A (en) * 2020-07-21 2020-11-24 沈阳安新自动化控制有限公司 Central air conditioner fault prediction and diagnosis method adopting double intelligent algorithms

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Research on State Monitor Techniques for Electronic System Based on HSMM;Li Wanling等;《2012 International Conference on Industrial Control and Electronics Engineering》;第1096-1098页 *
一种开关磁阻电机位置信号故障诊断与容错控制方法;胡荣光等;《电工技术学报》;第29卷(第7期);第104-113页 *
基于信号稀疏表示的滚动轴承微弱故障特征提取方法研究;张菀;《中国博士论文全文数据库》;工程科技Ⅱ辑 C029-16 *
基于动态HMM和离散SPC的空调机组的故障诊断;严颖等;《控制工程》;第29卷(第9期);第1707-1712页 *
基于概率密度函数的时滞依赖故障检测与诊断;陈龙等;《计算机工程与应用》;第54卷(第6期);第257-263页 *
数控机床的多性能与多故障关联矩阵模型研究;潘立峰;《中国优秀硕士学位论文全文数据库》;工程科技Ⅰ辑 B022-480 *

Also Published As

Publication number Publication date
CN113158546A (en) 2021-07-23

Similar Documents

Publication Publication Date Title
Zou et al. Towards optimal control of air handling units using deep reinforcement learning and recurrent neural network
Wang et al. Online model-based fault detection and diagnosis strategy for VAV air handling units
Nassif et al. Self-tuning dynamic models of HVAC system components
Jiang et al. Applying grey forecasting to predicting the operating energy performance of air cooled water chillers
CN114484731B (en) Central air conditioner fault diagnosis method and device based on stacking fusion algorithm
Yan et al. Fault prognosis of HVAC air handling unit and its components using hidden-semi Markov model and statistical process control
Li et al. Multi-objective optimization of HVAC system using NSPSO and Kriging algorithms—A case study
Yan et al. Fault diagnosis of components and sensors in HVAC air handling systems with new types of faults
Sun et al. Flow measurement uncertainty quantification for building central cooling systems with multiple water-cooled chillers using a Bayesian approach
CN113587172A (en) Water supply temperature delay time prediction method and device and electronic equipment
Homod et al. An innovative clustering technique to generate hybrid modeling of cooling coils for energy analysis: A case study for control performance in HVAC systems
CN113158546B (en) Fault prediction method of air conditioning unit equipment based on HSMM
Kim et al. Development of flow rate and equipment simulation model for commercial building HVAC&R system by data-driven method
CN113188818B (en) Air conditioning unit fault diagnosis method based on DHMM and CPDSPC
US11580281B2 (en) System and method for designing heating, ventilating, and air-conditioning (HVAC) systems
Nassif Modeling and Optimization of HVAC Systems Using Artificial Intelligence Approaches.
TW202204832A (en) Air conditioning system control method configured to control an air-conditioning setting temperature of an air-conditioning system and implemented by an air-conditioning control device
Adetola et al. Model predictive control and fault detection and diagnostics of a building heating, ventilation, and air conditioning system
Sun et al. Research on fault detection method for air handling units system
JP4434905B2 (en) Process management device
Huang et al. Enhancing the reliability of chiller control using fused measurement of building cooling load
Heo et al. Self-Training of a Fault-Free Model for Residential Air Conditioner Fault Detection and Diagnostics
US11162701B2 (en) Controlling HVAC system by inversing airflow dynamics
Wen et al. AHU AFDD
CN113029533B (en) Cooling coil fault diagnosis method based on hidden Markov model and discrete P control diagram

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant