Development of control quality factor for HVAC control loop performance assessment I—Methodology (ASHRAE RP-1587)

Abstract

This article is the second paper from the ASHRAE research project RP-1587, focusing on the methodology of obtaining a control quality factor (CQF). This article presents a development of two CQFs for assessing the heating, ventilation, and air-conditioning (HVAC) control loop performance. Both CQFs are able to detect whether the control loop is able to maintain the setpoint and identify the loop’s ability to handle disturbances. In addition, the reversal behaviors are assessed as well. The first CQF (CQF-Harris) is proposed based on the normalized Harris index using the recursive least squares method. This recursive least squares method is selected because of its computational efficiency compared with the maximum likelihood estimation method. The second CQF (CQF-EWMA) is based on the exponentially weighted moving average of the error ratio. The assessment scale of excellent, good, fair, bad, and failed, which indicates the quality of the HVAC control loops, is established as well. The sensitivity analysis for both CQFs is also conducted, and it provides insights on choosing the appropriate parameters to compute such CQFs. Such parameters include the sampling frequency, the length of the moving window, and the variance of the unmeasured disturbance. The field evaluations and tests of the proposed CQFs for simulated control loops and real control loops can be found in the companion paper with the title “Development of Control Quality Factor for HVAC Control Loop Performance Assessment—III: Field Testing and Results (ASHRAE RP-1587).”

Acknowledgments

Part of this work was supported by ASHRAE through RP-1587 “Control Loop Performance Assessment.” This work emerged from the Annex 60 project, an international project conducted under the umbrella of the International Energy Agency (IEA) within the Energy in Buildings and Communities (EBC) Programme. Annex 60 developed and demonstrated new generation of computational tools for building and community energy systems based on Modelica, Functional Mockup Interface, and BIM standards.

Appendix A. ARMA model for HVAC control loop output

For the ARMA fitting, the order of AR and MA modeling is very important. Numeral approaches have been proposed to select the appropriate order to minimize the fitting error. The widely used approaches are the autocorrelation plot (Box and Jenkins 1976), partial autocorrelation plot (Box and Jenkins 1976), Akaike information criterion with correction (AICc) (Hurvich and Tsai 1989), and Bayesian information criterion (BIC) (Schwarz 1978). However, they all have some drawbacks. This is still an active research area. There are some new progress lately for the order selection for ARMA models, such as the focused information criterion (Claeskens and Hjort 2003; Rohan and Ramanathan 2011). In our research, we select the order of AR to be 2, and the order of MA to be 1, with acceptable accuracy.

After we decide the correct ARMA order, the next step is to get the impulse response coefficient. Here we derive the equations on how to convert the ARMA model into the MA model, from which we can get the impulse response coefficients. The lag operator (denoted as L) is usually applied to simplified the converting. The lag operator is used to move the process output one time step back. Multiple (k) use of L will move the process output back multiple (k) times. The details are shown next: (22) $$\mathit{L}{\mathit{y}}_{\mathit{t}}={\mathit{y}}_{\mathit{t}-\mathbf{1}}$$(22) (23) $${\mathit{L}}^{\mathit{k}}{\mathit{y}}_{\mathit{t}}={\mathit{y}}_{\mathit{t}-\mathit{k}}$$(23)

The lag operator is also applicable to the addition, such as (24) $$\mathit{L}\left({\mathit{y}}_{\mathit{t}}+{\mathit{x}}_{\mathit{t}}\right)={\mathit{y}}_{\mathit{t}-\mathbf{1}}+{\mathit{x}}_{\mathit{t}-\mathbf{1}}$$(24)

For the ARMA model equation, (25) $${\mathit{y}}_{\mathit{t}}={\varphi }_{\mathbf{1}}{\mathit{y}}_{\mathit{t}-\mathbf{1}}+{\varphi }_{\mathbf{2}}{\mathit{y}}_{\mathit{t}-\mathbf{2}}+\dots +{\varphi }_{\mathit{p}}{\mathit{y}}_{\mathit{t}-\mathit{p}}+{\epsilon }_{\mathit{t}}+{\mathit{\theta }}_{\mathbf{1}}{\epsilon }_{\mathit{t}-\mathbf{1}}+{\mathit{\theta }}_{\mathbf{2}}{\epsilon }_{\mathit{t}-\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\epsilon }_{\mathit{t}-\mathit{q}}$$(25)

We rewrite the ARMA equation of the process outputs (Equation 25(25) $${\mathit{y}}_{\mathit{t}}={\varphi }_{\mathbf{1}}{\mathit{y}}_{\mathit{t}-\mathbf{1}}+{\varphi }_{\mathbf{2}}{\mathit{y}}_{\mathit{t}-\mathbf{2}}+\dots +{\varphi }_{\mathit{p}}{\mathit{y}}_{\mathit{t}-\mathit{p}}+{\epsilon }_{\mathit{t}}+{\mathit{\theta }}_{\mathbf{1}}{\epsilon }_{\mathit{t}-\mathbf{1}}+{\mathit{\theta }}_{\mathbf{2}}{\epsilon }_{\mathit{t}-\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\epsilon }_{\mathit{t}-\mathit{q}}$$(25) ) into (26) $${\mathit{y}}_{\mathit{t}}-{\varphi }_{\mathbf{1}}{\mathit{y}}_{\mathit{t}-\mathbf{1}}-{\varphi }_{\mathbf{2}}{\mathit{y}}_{\mathit{t}-\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{y}}_{\mathit{t}-\mathit{p}}={\epsilon }_{\mathit{t}}+{\mathit{\theta }}_{\mathbf{1}}{\epsilon }_{\mathit{t}-\mathbf{1}}+{\mathit{\theta }}_{\mathbf{2}}{\epsilon }_{\mathit{t}-\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\epsilon }_{\mathit{t}-\mathit{q}}$$(26)

Applying the lag operator, we get the new equation: (27) $${\mathit{y}}_{\mathit{t}}\left(\mathbf{1}-{\varphi }_{\mathbf{1}}\mathit{L}-{\varphi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{L}}^{\mathit{p}}\right)={\epsilon }_{\mathit{t}}\left(\mathbf{1}+{\mathit{\theta }}_{\mathbf{1}}\mathit{L}+{\mathit{\theta }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\mathit{L}}^{\mathit{q}}\right)$$(27)

Further derivation can be written as: (28) $${\mathit{y}}_{\mathit{t}}={\left(\mathbf{1}-{\varphi }_{\mathbf{1}}\mathit{L}-{\varphi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{L}}^{\mathit{p}}\right)}^{-\mathbf{1}}\left(\mathbf{1}+{\mathit{\theta }}_{\mathbf{1}}\mathit{L}+{\mathit{\theta }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\mathit{L}}^{\mathit{q}}\right){\epsilon }_{\mathit{t}}$$(28)

To simplify the preceding equation the process outputs can be written as a MA model as follows: (29) $${\mathit{y}}_{\mathit{t}}={\mathbf{\Phi }}^{-\mathbf{1}}\left(\mathit{L}\right)\mathbf{\Theta }\left(\mathit{L}\right){\epsilon }_{\mathit{t}}$$(29) where (30) $$\mathbf{\Phi }\left(\mathit{L}\right)=\mathbf{1}-{\varphi }_{\mathbf{1}}\mathit{L}-{\varphi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{L}}^{\mathit{p}}$$(30) (31) $$\mathbf{\Theta }\left(\mathit{L}\right)=\mathbf{1}+{\mathit{\theta }}_{\mathbf{1}}\mathit{L}+{\mathit{\theta }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\mathit{L}}^{\mathit{q}}$$(31)

The ARMA model is converted to the new MA model as (32) $${\mathit{y}}_{\mathit{t}}=\mathbf{\Psi }\left(\mathit{L}\right){\epsilon }_{\mathit{t}}$$(32) in which the new MA coefficients are written as (33) $$\mathbf{\Psi }\left(\mathit{L}\right)={\mathbf{\Phi }}^{-\mathbf{1}}\left(\mathit{L}\right)\mathbf{\Theta }\left(\mathit{L}\right)$$(33)

Expand the new MA coefficients with lag operators as (34) $$\mathbf{\Psi }\left(\mathit{L}\right)={\mathbf{\phi }}_{\mathbf{0}}+{\mathbf{\phi }}_{\mathbf{1}}\mathit{L}+{\mathbf{\phi }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathbf{\phi }}_{\mathit{s}}{\mathit{L}}^{\mathit{s}}$$(34)

Next, we obtain the new MA coefficients from the original AR coefficients and MA coefficients by (35) $$\mathbf{\Phi }\left(\mathit{L}\right)\mathbf{\Psi }\left(\mathit{L}\right)=\mathbf{\Theta }\left(\mathit{L}\right)$$(35) (36) $$\mathbf{\Theta }\left(\mathit{L}\right)=\mathbf{\Phi }\left(\mathit{L}\right)\mathbf{\Psi }\left(\mathit{L}\right)$$(36)

Expand the coefficients: (37) $$\mathbf{1}+{\mathit{\theta }}_{\mathbf{1}}\mathit{L}+{\mathit{\theta }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\mathit{L}}^{\mathit{q}}=[{\phi }_{\mathbf{0}}+{\phi }_{\mathbf{1}}\mathit{L}+{\phi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathbf{\phi }}_{\mathit{s}}{\mathit{L}}^{\mathit{s}}][\mathbf{1}-{\varphi }_{\mathbf{1}}\mathit{L}-{\varphi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{L}}^{\mathit{p}}]\mathrm{}$$(37)

Equation 37(37) $$\mathbf{1}+{\mathit{\theta }}_{\mathbf{1}}\mathit{L}+{\mathit{\theta }}_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathit{\theta }}_{\mathit{q}}{\mathit{L}}^{\mathit{q}}=[{\phi }_{\mathbf{0}}+{\phi }_{\mathbf{1}}\mathit{L}+{\phi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}+\dots +{\mathbf{\phi }}_{\mathit{s}}{\mathit{L}}^{\mathit{s}}][\mathbf{1}-{\varphi }_{\mathbf{1}}\mathit{L}-{\varphi }_{\mathbf{2}}{\mathit{L}}^{\mathbf{2}}-\dots -{\varphi }_{\mathit{p}}{\mathit{L}}^{\mathit{p}}]\mathrm{}$$(37) can be solved for the ARMA (2, 2) model without losing generality, and we can get: (38) $$\mathbf{1}={\phi }_{\mathbf{0}}$$(38) (39) $${\mathit{\theta }}_{\mathbf{1}}={\phi }_{\mathbf{1}}-{\varphi }_{\mathbf{1}}{\phi }_{\mathbf{0}}$$(39) (40) $${\mathit{\theta }}_{\mathbf{2}}={-\varphi }_{\mathbf{2}}{\phi }_{\mathbf{0}}-{\varphi }_{\mathbf{1}}{\phi }_{\mathbf{1}}+{\phi }_{\mathbf{2}}$$(40)

Thus, the impulse response coefficients are obtained as: (41) $${\phi }_{\mathbf{0}}=\mathbf{1}$$(41) (42) $${\phi }_{\mathbf{1}}={\mathit{\theta }}_{\mathbf{1}}+{\varphi }_{\mathbf{1}}{\phi }_{\mathbf{0}}$$(42) (43) $${\phi }_{\mathbf{2}}={\mathit{\theta }}_{\mathbf{2}}+{\varphi }_{\mathbf{2}}{\phi }_{\mathbf{0}}+{\varphi }_{\mathbf{1}}{\phi }_{\mathbf{1}}$$(43)

In all, the calculated impulse response coefficients are used for Harris index calculation.

Appendix B. Recursive least square method

The RLS is similar to the MLE method. Instead, it takes another approach by minimizing the cost function. This part presents the principles of the recursive least square method for ARMA fitting. We use ARMA(2, 2) as an example to show the whole process.

For an ARMA(2, 2) model, we have the following equation: (44) $$y\left(t\right)=a_{\mathbf{1}}y\left(t-\mathbf{1}\right)+a_{\mathbf{2}}y\left(t-\mathbf{2}\right)+b_{\mathbf{1}}w\left(t-\mathbf{1}\right)+b_{\mathbf{2}}w\left(t-\mathbf{2}\right)+w\left(t\right)$$(44)

Rewrite Equation 44(44) $$y\left(t\right)=a_{\mathbf{1}}y\left(t-\mathbf{1}\right)+a_{\mathbf{2}}y\left(t-\mathbf{2}\right)+b_{\mathbf{1}}w\left(t-\mathbf{1}\right)+b_{\mathbf{2}}w\left(t-\mathbf{2}\right)+w\left(t\right)$$(44) into the vector format: (45) $$y\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}y(t-\mathbf{1})& y(t-\mathbf{2})\end{array}& \begin{array}{cc}w(t-\mathbf{1})& w(t-\mathbf{2})\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{a}_{\mathbf{1}}\\ a_{\mathbf{2}}\end{array}\\ \begin{array}{c}{b}_{\mathbf{1}}\\ b_{\mathbf{2}}\end{array}\end{array}\right]+w\left(t\right)$$(45)

We get the simplified vector form (46) $$y\left(t\right)=h\left(t\right)\theta \left(t\right)+w\left(t\right)$$(46) where (47) $$h\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}y(t-\mathbf{1})& y(t-\mathbf{2})\end{array}& \begin{array}{cc}w(t-\mathbf{1})& w(t-\mathbf{2})\end{array}\end{array}\right]$$(47) (48) $$\theta \left(t\right)=\left[\begin{array}{c}\begin{array}{c}{a}_{\mathbf{1}}\\ a_{\mathbf{2}}\end{array}\\ \begin{array}{c}{b}_{\mathbf{1}}\\ b_{\mathbf{2}}\end{array}\end{array}\right]$$(48) and w(t) is the unmeasured disturbance term at current time step.

For a specified moving window with a sample of n data, we get the vector format: (49) $$Y\left(t\right)=H\left(t\right)\theta \left(t\right)+W\left(t\right)$$(49) (50) $$Y\left(t\right)=\left[\begin{array}{c}\begin{array}{c}y\left(t\right)\\ y(t-\mathbf{1})\end{array}\\ \begin{array}{c}\dots \\ y(t-n)\end{array}\end{array}\right]$$(50) (51) $$W\left(t\right)=\left[\begin{array}{c}\begin{array}{c}w\left(t\right)\\ w(\mathrm{t}-\mathbf{1})\end{array}\\ \begin{array}{c}\dots \\ w(t-n)\end{array}\end{array}\right]$$(51) (52) $$H\left(t\right)=\left[\begin{array}{c}\begin{array}{c}h\left(t\right)\\ h(t-\mathbf{1})\end{array}\\ \begin{array}{c}\dots \\ h(t-n+\mathbf{1})\end{array}\end{array}\right]$$(52) (53) $$H\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}y(t-\mathbf{1})\\ y(t-\mathbf{2})\end{array}& \begin{array}{c}y(t-\mathbf{2})\\ y(t-3)\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}w(t-\mathbf{1})\\ w(t-\mathbf{2})\end{array}& \begin{array}{c}w(t-\mathbf{2})\\ w(t-3)\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{c}\dots \\ y(t-n+\mathbf{1})\end{array}& \begin{array}{c}\dots \\ y(t-n)\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\dots \\ w(t-n+\mathbf{1})\end{array}& \begin{array}{c}\dots \\ w(t-n)\end{array}\end{array}\end{array}\right]$$(53)

According to the least square principles, we need to minimize the function: (54) $$J\left(\theta \right)=\frac{\mathbf{1}}{\mathbf{2}}\sum _{i=\mathbf{1}}^n{\left[y\left(i\right)-h\left(i\right)\theta \left(i\right)\right]}^{\mathbf{2}}$$(54)

The recursive steps updating the parameters are given as (LJung 1999): (55) $$\widehat{\theta }\left(t\right)=\mathrm{}\widehat{\theta }\left(t-\mathbf{1}\right)+P\left(t\right)\widehat{h}\left(t\right)\left[y\left(t\right)-\widehat{h}\left(t\right)\widehat{\theta }\left(t-\mathbf{1}\right)\right]$$(55) (56) $$P\left(t\right)=P\left(t-\mathbf{1}\right)-\frac{P(t-\mathbf{1})\widehat{h}\left(t\right){\widehat{h}}^T\left(t\right)P\left(t-\mathbf{1}\right)}{\mathbf{1}+{\widehat{h}}^T\left(t\right)P\left(t-\mathbf{1}\right)\widehat{h}\left(t\right)}$$(56) (57) $$\widehat{w}\left(t\right)=y\left(t\right)-\widehat{h}\left(t\right)\widehat{\theta }\left(t\right)$$(57) (58) $$\widehat{h}\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}y(t-\mathbf{1})& y(t-\mathbf{2})\end{array}& \begin{array}{cc}\widehat{w}\left(t-\mathbf{1}\right)& \widehat{w}\left(t-\mathbf{2}\right)\end{array}\end{array}\right]$$(58)

The initial conditions are given as: (59) $$P\left(\mathbf{0}\right)=P_{\mathbf{0}}I$$(59) (60) $$\theta \left(\mathbf{0}\right)=I_n/P_{\mathbf{0}}$$(60) (61) $$P_{\mathbf{0}}={10}^6$$(61)

The initial condition of y(1) is given as: (62) $$y\left(\mathbf{1}\right)=a_{\mathbf{1}}y\left(\mathbf{0}\right)+a_{\mathbf{2}}y\left(-\mathbf{1}\right)+b_{\mathbf{1}}w\left(\mathbf{0}\right)+b_{\mathbf{2}}w\left(-\mathbf{1}\right)$$(62) (63) $$y\left(\mathbf{1}\right)=\left[\begin{array}{cc}\begin{array}{cc}y\left(\mathbf{0}\right)& y(-\mathbf{1})\end{array}& \begin{array}{cc}w\left(\mathbf{0}\right)& w\left(-\mathbf{1}\right)\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{a}_{\mathbf{1}}\\ a_{\mathbf{2}}\end{array}\\ \begin{array}{c}{b}_{\mathbf{1}}\\ b_{\mathbf{2}}\end{array}\end{array}\right]=h\left(\mathbf{0}\right)\theta \left(\mathbf{0}\right)$$(63)

In all, we obtain the AR and MA coefficients based on this RLS method. The AR and MA coefficients are used for the Harris index calculation.