Korkut Bekiroglu, Seshadhri Srinivasan, Ethan Png, Rong Su,, and Constantino Lagoa

Abstract—This paper presents an approach to recursively estimate the simplest linear model that approximates the time-varying local behaviors from imperfect (noisy and incomplete) meas-urements in the internet of things (IoT) based distributed decision-making problems. We first show that the problem of finding the lowest order model for a multi-input single-output system is a cardinality ? optimization problem, known to be NP-hard.To solve the problem a simpler approach is proposed which uses the recently developed atomic norm concept and the modified Frank-Wolfe (mFW) algorithm is introduced. Further, the paper computes the minimum data-rate required for computing the models with imperfect measurements. The proposed approach is illustrated on a building heating, ventilation, and air-conditioning(HVAC) control system that aims at optimizing energy consumption in commercial buildings using IoT devices in a distributed manner. The HVAC control application requires recursive thermal dynamical model updates due to frequently changing conditions and non-linear dynamics. We show that the method proposed in this paper can approximate such complex dynamics on single-board computers interfaced to sensors using unreliable communication channels. Real-time experiments on HVAC systems and simulation studies are used to illustrate the proposed method.

I. Introduction

A. Motivation

SCIENTISTS have recently put significant effort in decentralizing most of these centralized problems to improve the general performance of the systems [1]. Furthermore, in many distributed decision-making problems, the cost of computation may escalate because of the structure of the local agents’models. Moreover, complex models require exhaustive computations for control systems applications. Therefore, implementable distributed control of complex processes requires the use of “simple models” that can i) provide good local approximations of the behavior; and ii) lead to control algorithms that are computationally efficient and easily implementable.Such problems are widespread in the internet of things (IoT),due to limited computing resources, imprecise sensors, and data imperfections resulting from unreliable networks [2].Also, adaptability is an essential feature in many emerging applications, including IoT, and therefore, recursive identification becomes pivotal for such systems. The problem becomes interesting in the building heating, ventilation, and air-conditioning systems due to the multiple thermal zones in a building that have complex non-linear behaviors and subjected to frequent disturbances [3], [4].

B. Literature Review

Data imperfection can arise in many applications. For instance, incomplete measurements arise in many applications such as fault detection in building HVAC systems [5]. In the system identification literature, subspace identification methods are widely used in both linear [6] and non-linear system identification [7]. However, these techniques cannot handle sparsely sampled data, and furthermore, they are not amenable to recursive implementation. Therefore, considerable effort has been recently devoted to handling missing measurements in system identification problems while simultaneously enforcing system order [8], [9]. Most of these methods use a rank minimization (RM) method, a non-convex NP-hard problem that is mostly relaxed by using the concept of nuclear norm [9], [10]. In this scenario, the missing measurement is a decision variable of the optimization problem which compounds the computation complexity further.Recently, a method to identify multi-input and single-output(MISO) models with missing output data was proposed in[11]. However, the recursive identification was not considered in the presence of missing measurements. In addition to linear system identification, non-linear system identification literature also has many mature methods [7] to approximate the non-linear behavior of systems. Even though these methods have recursive model identification, according to our best knowledge literature still needs local low order linear model approximation of non-linear systems from the incomplete dataset. As it is mentioned in the motivation Section I-A, having a local low order linear model approximation is an essential step for distributed control applications such as model predictive control (MPC) to diminish the complexity.

C. Contributions

In this paper, we present an algorithm which recursively computes simple linear models that approximate complex non-linear and time-varying local behaviors in the presence of imperfect (noisy and incomplete) measurements. A review of the literature reveals the lack of techniques for computing lower-order system models recursively from input-output data with missing measurements, a common occurrence in many distributed systems. In contrast to existing results, this paper presents an approach that directly identifies the lowest order MISO model recursively in the presence of imperfect measurements. We should note that the MIMO version of the proposed work requires further research. The proposed method only assumes that the system is stable and persistently excited, unlike existing approaches that require additional information on model order, accuracy, and continuous data.The main contributions of the investigation can be summarized as:

1) A method to recursively identify complex non-linear and time-varying behaviors using simple (low order) MISO linear models for distributed decision-making problems.

2) Compute the minimum data-rate required for identifying such lower order models from imperfect measurements.

3) Illustration of the proposed method in a pilot IoT based building management system which shows the efficiency and performance of the proposed approach.

D. Organization

The problem statement is presented in Section II. The relaxation for the cardinality optimization problem and a modified Frank-Wolfe algorithm to solve the relaxed problem are presented in Section III. Section IV presents the theory for the minimum data requirement for system identification.Illustrative examples are provided in Section V. Conclusions and further directions of this investigation are presented in Section VI.

Notations:

II. Problem Statement

Consider the distributed decision-making problem with each controlled entities having the MISO dynamical model

where F represents a non-linear time-varying function, y(tk),u1(tk),u2(tk),...,uN(tk) , and ε(t) denotes the sensor measurements (output), inputs, and the measurement noise, respectively, andny,nuare the time lags for the system output and inputs. Our idea is to approximate the non-linear and time-varying model in (1) of each controlled entities with the simplest linear model which captures the essential dynamics with imperfect sensor measurements but known inputs1The imperfect measurements can result from a variety of reasons which include hardware aspects such as sensor failures, noise, faults, etc., or network induced imperfections such as latencies, packet dropout, quantization, and others.. The following definitions are given to formulating the proposed problem.

Definition 1 (Linear Local Measurement Model):Considering y(tk) to be the measurements obtained with commensurate sampling instantstk:t0

where * denotes convolution,mis the number of measurements, girepresents the impulse response from inputito output, uirepresents the input signal ofith input, and y(tk) is output signal at timetkcorrupted by some bounded noise ε.

Definition 2 (Stable Dynamic Region of the Model):The dynamics of the simple linear model are characterized by a set of poles ?i∈Z+,pi∈Dρ, inside a closed centered disc around the origin with ρ ≤1. The set of poles are deemed as well behaved, i.e.,pi∈Dρ={pi∈C orpi∈R:|pi|≤ρ}2The set of poles can be defined based on priori such as rise time, overshoot,etc. Only the stable poles are chosen for this paper, reasonable one for HVAC systems, but can include unstable poles..

Definition 3 (Data Loss Operator):Data-loss is modelled by a binary variable σtk∈{0,1}, with σtk=1, if the packet is delivered and 0, otherwise. The network communication formpackets is modelled as am×mmatrix Lm=diag[σt0...σtm],and the channel is assumed to be perfect, iftrace(Lm)=m+1,and imperfect otherwise. The data loss operator models the constraints due to data imperfections in the model identification problem. It is a LmBmatrix multiplication operator which models the data loss. The operator takes as input the tuple {Lm,B} and provides an output that replaces the rows ofBwith zeros for the rows with all zeros in the matrix Lmand leaves the other rows unaltered.

Without loss of generality, the measurement model in (2)can be extended to model the IoT data constraint by introducing an IoT data loss operator Lm, i.e.,

Remark 1:The linear approximation of the non-linear and time-varying systems is ill-defined. Given a slow time-varying non-linear system, one can locally (in-state and time)approximate its behavior by an LTI system. Since there are many ways to approximate this behavior, we choose an approximation that provides a trade-off between complexity and error.

Finally, we tried to address the following system identification problem in Problem 1 under the following Assumptions 1–3.

Problem 1:Suppose there exists a lowest order modelgithat approximates the non-linear and time-varying local behaviour of a system. The problem is to computean estimate of the gihaving poles in the stable dynamic regionDρ(see, Definition 2) from imperfect measurementsunder the given assumptions.

Assumption 1:The lowest order transfer functionhas no repeated poles. This assumption does not limit the capability of the proposed method since any repeated poles can be approximated with non-repeated poles, but with an arbitrary small perturbation [12].

Assumption 2:Initial conditions are assumed to be zero to simplify the analysis of the proposed method. Non-zero initial conditions can also be addressed with a small modification of the optimization problem3Discrete linear systems with distinct poles have two facts: i) the same poles appear in the impulse response and the initial condition response; ii) the number of distinct poles used in these two is equal to the order of the system.Hence, sparsifying the impulse response of the system (i.e., identifying the system of lowest order) is equivalent to minimizing the number of unique poles that are used in the total response, i.e., the response that includes both the response to initial conditions and the zero state response to the input..

Assumption 3:The measurement noise in (2) is bounded,i.e., ∥ ε∥2≤Ξ.

In what follows, we show that the problem of finding the lowest order linear model recursively with imperfect measurements that captures complex non-linear behavior of each of the controlled entities is a ?0optimization model. We propose an elegant way and an efficient algorithm to solve the complex problem. Furthermore, we also investigate the data that needs to be transmitted for approximating the non-linear behaviours.

III. Recursive System Identification Methodology With Imperfect Data

In this section, the optimization problem of the minimum order system identification of a MISO system is presented.

A. Minimum Order Model Identification Problem

The impulse response of a linear system giin (2) can be expressed as a linear combination of a set of first-order systems’ impulse responses, denoted byAtoms, by using partial fraction expansion [13]. To obtain the minimum order transfer function, the cardinality ofAtoms(number of these first-order systems) has to be minimized. However, one should note that even though each transfer function of inputito output shares the same poles, they could have different zeros. This requirement can be addressed using the following approach

where ?0norm minimization of a vector c in (4a) is analogous to minimizing cardinality (c ar dinality{c:}) wherecpis thepth element of vector c. The constraint in (4b) enforces the accuracy of the model for a bounded noise for real output data y(note that we assume complete data in this step and incomplete data will be address in (13)) and (4c) enforces a block sparsity ofwhich means the sparsity pattern of each input vector are same. In other words, it is ensured that each transfer function from output to inputiuse same poles. Note thatrepresents a constant for each first order system for inputi.There are four challenges in the optimization problem in (4a),they are: i) the cardinality minimization problem is a NP-hard problem [14]; ii) the computation complexity arising frombeing complex; iii) the search-space (Dρ) contains infinite number of poles; and iv) the need to perform recursive identification compounds the complexity further.

B. Convex Relaxation

To relax the optimization problem in (4a), this paper utilizes an atomic norm-based approach. The complete definition and properties of the atomic norm can be found in [15]. From the definition ofAtoms, if one chooses simple objects asAtomsin an atomic set and tries to force using few of them in expressing a certain object, this naturally promotes a sparse representation of the object at hand. The proposed method is developed by exploiting this property. In our analysis, we choose the impulse response as theAtomsas given in Appendix-A.

The noisy measurement model y(tk) of the system (2) can now be represented using the following matrix and vector notations as

where Tui∈Rm×mis the Toeplitz matrix ofith input and giis the impulse response for the inputuito output y ∈Rm. The impulse response of the transfer function giin (5) is approximated as

where Z?1denotes the inversez-transform. The sparsity pattern of eachin (7) withbeing theith element of vector cais identical to enforce the same pole in each of the transfer function,gi

Rewriting (5), we have

with

and

One can realize that each giutilizes the same atom. Using the reformulation in (8), the non-convex optimization problem in (4a) can be relaxed using the following atomic norm

Having defined the atomic norm that encourages block sparsity for the MISO system transfer function, the following convex constrained optimization problem (12), a convex relaxation problem of (4a), can be solved. This relaxes the requirement of solving an NP-hard problem.

where τ is a bound on the convex surrogate of block sparsity of vectors c in (4a). The IoT induced constraints are brought in by introducing the data loss operator Lmas

With this trivial modification, the optimization model for obtaining the minimum order transfer function for the distributed decision-making system with imperfect data is obtained. Though the non-convex ?0optimization problem is now relaxed to a convex problem, the presence of infinitely many poles in the setDρincreases the search space significantly. The presence of infinitely many poles is solved in two ways: griding the region of the well-behaved poles [16]and a randomized Frank-Wolfe algorithm (rFWA) [11]. While griding is a computationally inefficient method to search, the rFWA is not suitable for recursive computations. Therefore,this investigation introduces a heuristic step that enhances the computational speed and makes the algorithm suitable for recursive system identification. We call this modified version of rFWA as a modified Frank Wolfe algorithm (mFWA).

C. Solving Convex Optimization by Using mFWA

In this section, we introduce a method to solve the convex problem in (13) that is a special case of the atomic norm constraint problem presented in [17];

wheref(g) is a convex and smooth function. In our case, the functionf(·) is thei.e., the 2-norm of the error and the definition of ∥ g∥Aare given in (11).

The various steps involved in the mFWA are shown in Algorithm 1. The inputs to the algorithm are the initial modelthe compact pole location setDt, bounds on the iteration error ?, the convergence criteria γ , τ bound on the atomic norm of the impulse response and data, and the pole locations of initial model Dt,Pt. The convergence criteria, iteration bounds, and iteration counterkare first initialized. The following computations are performed until the convergence criteria and iteration bounds are met.

i) The polesare selected by uniformly samplingDt, a region having stable poles (Step 4).

iii) The gradient vector of the objective function in (13) is calculated as

This gradient vector is partitioned for each impulse responseiasStep 6).

iv) Different from rFWA, we modified this step in our algorithm for MISO systems. The minimum over the set of the inner product of the gradient and the + and ? signs of correspondingAtomsis computed for each inputi. This gives the descent direction for the optimization problem (Step 6).

Remark 2 (Promoting Block Sparsity):The gradient vector for each impulse responseis treated with an atom±ain atomic setThis enforces each impulse response to use the same pole. However, because of the structure of the atomic setAin (22), the sign combination of this chosen impulse response is tested to find the descent direction for each impulse response.

v) The optimal step-size αkfor each iterationkwith givenτ andfor the best atomakwith α ∈[0,1] is computed as (Step 7)where

where

where

and

vi) Then this step is different than the classical rFWA for SISO systems since this step is added to promote the recursive system identification. In here, if thethen the selected atom is used to update the impulse response. Therefore, the pole corresponding to the atom, i.e., (pk) is appended to the existing pole set Ptfor the recursive computationt+1 (Step 8).

vii) The following heuristic is also added to rFWA to improve the quality of the algorithm for the MISO system identification. Although the bound τ on the atomic norm ofenforces the block sparsity of vectorsin (7), magnitude of each impulse response gishould be adapted in each iteration with an additional constant. Therefore, the scaling factoris calculated in each step to scale the magnitude of impulse response from inputito output (Step 10). After calculating each impulse response, the impulse responseis scaled with(Step 9)

viii) Update the impulse response of each inputiasand compute the matrix(Step 11).

ix) The loop is repeated until convergence. The output of the algorithm is the

Algorit hm 1: mFWA

Algorit hm 1: mFWA

For the rFWA, it has been already proven that rFWA converges in expectation [13]. In addition to convergence in expectation, the rFWA for SISO case convergences with a fixed rate. Since the randomization step in rFWA and mFWA is the same (Steps 4 and 6), we have the same rate of convergence for our mFWA. For the completeness, we highlighted Lemma 1 and the Theorem 2 given in Appendix-B for our recursive identification algorithm.

In what follows, we will illustrate how the mFWA proposed in this investigation can be used for recursive system identification.

D. Recursive Minimum Order System Identification Methodology

In this section, we introduce one of our main contributions.The proposed recursive system identification methodology is shown in Fig. 1. The method works in two modes: off-line and on-line. In the off-line mode, the data is used to learn a minimum order nominal modelfrom a set of data D0by solving (13). In the off-line mode, the bound on the approximation error, an input of the Algorithm 2, can be also learned by analyzing the fitting error. This error is the approximation errors of linear estimation of non-linear and time-varying behavior. In this way, the user can have an initial guess of τ in (14). Then this model (impulse response)obtained in the off-line mode and its poles are inputs to the on-line mode. The heart of the model adaptation algorithm is the mFWA which takes the nominal model from the previous step, the current set of data D1which also has missing measurements, and the pole locations of this nominal model to obtain the lowest order model. The procedure repeats when a new set of data is available. In this recursive system identification algorithm hot-starts the next iteration using the poles of the current model.

Fig. 1. Recursive identification methodology.

While in the off-line mode, the poles are randomly chosen inDρ(see, Definition 2), in the on-line mode, the poles are randomly chosen in the norm ball around each pole of the nominal model computed in the previous step. The norm-ball is defined as follows.

Definition 4 (Norm Ball Around a Pole p):Given the polep, a closed centered norm ball around the polepwith radiusβ is the set of poles?i∈Z+,pi∈Dβpsuch that{pi∈C orpi∈R with ∥pi?p∥≤β}.

The concept of norm-ball enables searching the poles in the on-line mode around theDβpregion instead of the whole regionDρ, and this increases the speed of computation significantly as the most likely locations are checked by the mFWA first. Our main algorithmic contribution is given in Algorithm 2.

Algorithm 2: Recursive Method

Furthermore, in between each mode, we include a model reduction step because of the structure of the Frank-Wolfe algorithm. One can see in Step 9 in Algorithm 1, the poles in each step are appended to the current impulse response. This process directly affects the parsimony of the system in the recursive iteration. Therefore, during the recursive system identification, some of the poles in the previous step are eliminated using the following method.

Remark 3 (Model Order Reduction):Before initiating the new identification step by utilizing the previous nominal model with a new set of data (see, Algorithm 1), this previously identified impulse response is used to construct a Hankel matrix. Since the rank of the Hankel matrix is equal to the system order, the singular value decomposition (SVD) is calculated and the tails of these singular values are eliminated.Since this method is a well-known method, we refer readers to the following studies [18]–[20].

IV. Minimum Data Required for Linear System Identification

After presenting the recursive system identification algorithm in the previous section, in this section, we analyze how much data we need to identify anth order system. In other words, we investigate how much missing data can we handle to identify the true approximated system. Prediction error of annth order system from (8) iswhere c ∈Rnand y,are with suitable dimensions. The multiplication ofand the unknown vector c is the linear combination ofnimpulse responses. Then the following optimization computes the corresponding impulse responses in

Assuming the input u is a quasi-stationary signal, the uniqueness of c?is guaranteed as long as the matrixhas full column rank and can be determined as the solution of. However the data loss operator Lmin (13) reduces the rank of the matrix. From the definition of the persistent excitation condition [21], one can state that the rank of the correlation matrixis equal to the dimension of the c.

Using this fact, one can propose the following theorem.

Theorem 1 (Necessary Condition for Identifiability):Suppose there aremimperfect measurement packets and our aim is to identify a system whose highest order isn, wherem>>n. The necessary conditions for identifiability of the system with the imperfect measurement is given by

The proof of the theorem is given in Appendix-C.

V. Examples

This section presents two case-studies, a real-time experiment, and simulations that demonstrate how complex behaviors can be approximated as simple models by using the approaches presented in this paper for distributed decision making. Both the examples concern the building with variable air volume (VAV)-based HVAC control, a complex nonlinear system. Even though imperfect measurement can be seen in most of the IoT based implementations, HVAC control is one of the good cases to test the proposed work [22]. Note that estimating the energy model for building automation is an important problem.

A. Case-Study: IoT-Based HVAC Control

1) Pilot Description:The pilot considered is the VAV controlled HVAC system located at S1B1 building, NTU,Singapore. It has 85 zones and 6 AHUs sharing a centralized chiller [23]. The test-bed has AHU that are manufactured by York and have flow sensors from ABB. The fan power is 11 kW induction motor and the centralized chillers’ flow is controlled using a valve. The fan supplies air to individual zones through a common duct. The fan speed is controlled using a speed control drive and the power varies as a square of the fan speed. The test-bed has control systems supplied by Siemens A/s to handle the AHU control and Quantum systems for the chiller control. Currently, the test-bed has a centralized control architecture with each zone’s thermal behavior having a nonlinear dynamics (bi-linear) due to coupling between state and control inputs. We place IoT sensors to collect local information as explained later and there are issues of data-loss due to communication imperfections. In the experiment, we describe data collection, sensors, communication and other aspects with which data is collected. Currently, a centralized controller is used for allocating the mass-flow rate to each of the zones by controlling the opening of VAV. The air-flow in the zones are coupled through mass-flow rate constraints and fan capacity. Following [24] the zone thermal dynamics is given by

whereci,ca, andNzdenote the thermal capacitance of the zonei, specific heat capacity of air [kJkg–1°C–1] and the number of zones, respectively. Additionally,m˙i, Xc, Xi,Ri, and Xoadenote the mass-flow rate in zonei[kgs–1], temperature of the cool air from AHU [°C], temperature of the zonei[°C],thermal resistance of zonei[kW°C–1], and ambient temperature [°C], respectively. The termQ˙ models the heat load rate and is influenced by factors such as occupancy which are time-varying. Consequently, model adaptation is required for capturing influences of such time-varying factors.

2) Distributed Control:Current building automation systems are centralized structures and use conventional control techniques such as proportional-integral controller or thermostats which are not optimal. Recent research has proven that using the information on weather, occupancy, etc., to control mass-flow rate in individual zones can enhance energy efficiency. In this backdrop, the model predictive controller which uses the zone thermal models, coupling in the AHU and chiller circuits to compute the optimal mass-flow rate has emerged as a promising solution. Implementing MPC with centralized architecture is complex due to the non-linear and time-varying zone thermal dynamics. Distributed control wherein the optimal mass-flow rates are computed using single board computers interfaced with sensors that offer scalability and simplicity. Further, by using the IoT such boards can be networked to give a distributed implementation.However, the non-linear and time-varying zone thermal dynamics makes the realization of MPC on a single-board computer complex. Therefore, simple linear models that approximate complex thermal dynamics recursively are required for realizing distributed HVAC control.

3) IoT-Based Implementation:To recursively approximate the non-linear zone thermal dynamics with simple linear models, the linear model presented in [24] and the method presented in this paper were used. To model the zone thermal dynamics, we conducted experiments to record the VAV mass-flow rate in l s–1was measured and it was converted to kgs–1using a suitable conversion factor. In addition,influences of weather were measured using a web-service and the occupancy was measured using CO2sensors. The temperature and humidity were measured with DHT11 sensors interfaced with the Raspberry Pi board and the sensor measurements were transmitted to the controller board that implements the MPC and our model approximation approach proposed in the paper. We use the NRF2401L wireless module to transmit the measurement by constructing the sensor module as shown in Fig. 2. The input measurements from these sensors are given in Fig. 3. The sensor module in Fig. 2 is used to measure the mass-flow rate from which control input could be computed as

Fig. 2. Sensor module with wireless interface.

Fig. 3. Mass-flow rate and ambient temperature recorded during the experiment.

whereis the cooling energy rate,is the zone mass-flow rate,Cρis the specific heat capacity of air,TsandTzare the supply air temperature, and zone temperature, respectively.Since the data is transmitted to the controller module using a wireless network, there are packet-losses which make the identification challenge. Our experiment showed that about 20?25% of the data-loss occurred due to wireless communication and sensors. The need to have the lowest order model to reduce computations adds another dimension to the complexity. The parameters used in our experiments are τ=3,ρ=0.999, and β=0.15. The measurements arem˙i,Xoa, and

4) Experiment:The experiments were performed on 28th Feb. 2018 for a period of 200 minutes. The mass-flow rate and ambient temperature measured are shown in Fig. 3. To collect the data, we used Apache24 (a HTTP server), PHP7+HTML5(web), and MySQL (database). We chose these platforms as they are open-source that provides free to use software without licensing constraints. During the experiments, the temperature set-point of the thermostat controlling the VAV is varied from 23 to 27 °C. The model adaptation when the temperature settings changed is reported here. During the experiments the occupancy, heating load and other factors were allowed to be changed and only the set-point changes were made during specific instants. The first data-set (see,Table I(a)) shows the nominal impulse response models and initial condition response estimation computed by using the off-line imperfect data. The identified initial condition response, impulse response with cooling energy and ambient temperature as inputs and zone temperature as the output are also shown in the second and third column of Table I. In addition, the fourth column shows the actual response versus the one computed by the identified model (initial condition plus input responses that are the convolution of impulse responses and inputs). The performance of the proposed method for the first set of data is then evaluated with mean squared errorwhich is 1.59 and presented in the last column of Table I. The off-line model is the nominal model of the proposed approach and the models are adapted during the on-line mode.

During the on-line mode, imperfect measurements are used to update the model recursively. Four iterations of the model updates are shown in the Table I (b)?(e) with different percentages of packet-loss. The nominal model is updated by flipping the mass-flow rate fromhigh to low, and the 1st batch of the online data is recorded. The impulse responses with the cooling energy and ambient temperature as inputs and the actual measurements as well as the predicted measurements from the identified model are also shown. One can see that with the model adaptation step the modeling error, i.e., MSE reduces to 0.0047. The model obtained in the 2nd step is updated in the 3rd step and approximately 20% packet dropout is considered. One can see that the model accuracy increases due to recursive identification by about 20%. The4th set of experiments are conducted with 30% data-loss (the highest anticipated loss from the network) and the error increases to 0.0002. Still, the modeling error is low and indicates the efficiency of the approach to recursively approximate complex non-linear behaviors. The 5th batch of data is collected from the experiments with packet-loss of 20%. But due to the recursive model updates, the MSE is 0.0025over 20 samples. This shows that the accuracy is sensitive to packet losses. Since the model is simple (lowest order and linear), the MPC can be designed to reduce the energy consumption in a distributed manner. Note that in this study, we only identify the nominal model from noisy and incomplete measurements. The identified models were not implemented yet in the controller design stage since we only focused on the identification part of this study.

TABLE I(a) Nominal Model From Off-Line Data; (b) Updated Model Using Nominal Model in First Online Mode; (c) Online Adapted Model From the Previous Step; (d) Online Adapted Model From the Previous Step; and (e) Updated Model With the Fifth Data-Set ( (a) to (e) are Numbered Top Row to Bottom Row)

B. Simulation Case-Study

In the literature, prediction error method (PEM) and subspace identification methods have been widely studied for system identification. Similarly, two approaches or routines are used for handling imperfect data:MissdataandMerge.TheMissdataroutine linearly interpolates the missing measurements, and a full data-set is used for identification. On the other hand, the routineMergeremoves the NAN values and combines data-sets into a set of “multi-experiment data object” which is used for identification. While the proposed algorithm does not require any filling or data-processing approaches such asMissdataandMergefor handling imperfect data, other approaches in the literature are limited by this requirement.

To illustrate this aspect, the discrete-time version of the non-linear zone thermal dynamics in (18) sampling period δ is given by [24]

where ui(tk)=m˙i(Xi(tk)?Xc) is used to enforce linear assumption on the system. Consequently, the non-linear dynamics are approximated by a first-order model which can be used to implement sophisticated control actions. Without loss of generality, the approach can work for a system of any order and not restrained to first-order ones, the current example is used to simplify our analysis. The system model considered is given by [25] with ξ1=0.64,ξ2=?2.64 , and ξ3=0.1, where Xis the room temperature, u is the cooling energy input,wis the outside temperature, and ε is the noise. Given the zone thermal model and the sampling rate is 15 minutes, the model is excited with a PRBS input with a length of 3000 minutes.Furthermore, the temperature measure based on this input is contaminated with following white noiseN ～(0,0.25×max(XTrueMeasurment)). Then 30% of the temperature measurements are randomly eliminated to generate data-loss on the output. The proposed approach is used to identify the model with parameter τ=1.9, to provide a good trade-off between error and model order. The impulse response is computed aswith inputs u andw. The model identified by the proposed method with imperfect measurement is=0.6782 ,=–2.34, and

C. Benchmark Comparison

The model identified is compared with the PEM and subspace approach (n4sid) by usingMergetechnique, and a comparison of the impulse responses obtained is shown in Fig. 4.It is evident from the impulse response that the proposed method performs better than existing approaches even with imperfect measurements. A comparison of the model fitting is provided in Table II which also illustrates the utility of the proposed approach.

Similarly, the robustness of the approach to data-loss is tested by randomly exiting the system 500 times with the following conditions:

1) the sampling rate is 15 minutes;

2) the model is excited with a PRBS input with the length of 2250 minutes;

3) white noise isN(0,0.15×max(XTrueMeasurment))

4) when 50%, 60%, 70%, 80%, 90%, 100% of temperature measurements are available;

5) τ1=4.5 and τ2=0.2.

Fig. 4. Impulse response estimates.

TABLE II Proposed Approach vs PEM and Sub-Space Identification

TABLE III Average Goodness of Fit for Different Level of Data-Loss

The robustness is studied by defining the Goodness of Fit measure, given by and the performance of the algorithm is tabulated in Table III.One can observe that the quality of fitting is non-linear with respect to data-loss and the average performance is still within reasonable limits though the maximum is sensitive to data loss. These results illustrate the robustness of the proposed approach vis-á-vis, data-loss.

VI. Conclusion

This investigation proposed an approach to recursively estimate complex non-linear behaviors using simple linear models from imperfect measurements for distributed decisionmaking problems. The method focused on the Internet of Things based control of building heating, ventilation, and airconditioning system that performs distributed decision making with unreliable networks. We first showed that the problem of computing the lowest order model is a ?0optimization problem. The influence of multiple inputs and imperfect data made the problem even more challenging. An approach to solve this challenging problem was proposed using the atomic norm concept and a modified Frank-Wolfe algorithm. Further,the minimum amount of data required to identify the lowest order model was also studied to handle imperfect measurements. The proposed approach was illustrated in a building HVAC system both with experimental and simulation examples. Our results showed the ability of the method to compute the lower order linear model from imperfect measurements. Further, the method showed better performance than existing methods such as prediction error method and sub-space identification that require restrictive assumptions. Extensions to multi-input and multi-output systems are the future course of this investigation. Also, the proposed work can be integrated with building operation,management, or predictive control.

Appendix

A. Atomic Set

In the literature, to address the issue of complex termsin (4a)), instead of taking linear combination of first order system responses, second order ones have been used [26], leading to real-valued coefficients. To include this modification, we define an atomic set which is denoted bywhich takes set of poles, namedor(see,Definition 4 in Section III-D) as an input and generates theAtoms. The set consists of the first, second order responses,and the one to represent the poles in the unit circle as

with the scaling factors, αps, are used as a weighting for each atom to normalize the Hankel nuclear norm of it in the setA[16], [26]4Justification of Hankel norm normalization can be found in [16], [26].. This leads to

The atomic set is

wheremis the length of measurement vector y ∈Rm.

B. Convergence Lemma and Theorem

The converges in expectation lemma is

Lemma 1:Algorithm 1 converges in expectation and almost surely (a.s.) for a fixed number ofAtomssayNksatisfyingNk≥1 for allk, i.e., letf?be the optimum of problem (12) or(13), then

and

Also, the convergence rate is

Theorem 2:LetC2>0,L>0 , and 0

where g?is an optimal solution of Problem (13).

The proof of Lemma 1 and Theorem 2 are omitted since it is given in [13].

C. Proof of Theorem 1

Proof:The correlation matrix of the output from (13) is

From the definition of the convolution sum, we have

Considering that the system is persistently exited by the quasi-stationary input u, the following necessary condition needs to holdrank(Ru)≥n[21]. Clearly, for the rank deficient matrix Lm, for the system to be persistently excited, the rank conditions can be expanded

Moreover, it is known that ifthenand Assumption 2, we have

Therefore, the persistent excitation condition leads to

Necessary condition for identifiability ofnth order system is given byrank(Lm)≥n.

Remark 4:This necessary condition is sufficient for noisefree measurements or one can find a sufficient condition for minimum data requirement for noisy measurement if the noise frequency is known.