\section{Introduction}\label{sec:intro}
	A dynamic model of a building's temperature is necessary for model-based control of building HVAC  (Heating Ventilation and Air Conditioning) systems~\cite{Privara:2013}. Due to the complexity of thermal dynamics, system identification from data is considered advantageous and there has been much work on it; see
	~\cite{penman1990second,bracha:2002,wang2006parameter,YL_TM_PB_CDC:2012,liwen:2014,harb2016development,fux2014ekf,kimcaiaribra:2016,hu2016building} and references therein. A particular challenge for model identification is that temperature is affected by large, unknown disturbances, especially the cooling load induced by the occupants. The occupant induced load refers to the heat gain directly due to the occupants' body heat and indirectly from lights and other equipment they use. Many system
	identification methods ignore these disturbances \cite{YL_TM_PB_CDC:2012,liwen:2014,harb2016development}, or use a specialized test building to measure the occupant induced load~\cite{penman1990second,wang2006parameter}.
	Ignoring the disturbance can produce highly errorneous results~\cite{kimcaiaribra:2016}.

	In this paper we propose a method to estimate a dynamic model as well as a transformed version of the unknown disturbances from easily measurable input-output data. The proposed method, which we call SPDIR\ (Simultaneous Plant and Disturbance Identifier with Sparisification) is based on solving a $\ell_1$-regularized least-squares problem. The $\ell_1$ penalty encourages the identified transformed disturbance to be sparse~\cite{HastieElements:2009}. The motivation for this is that the disturbance, which consists mostly of internal load due to occupants, is often piecewise-constant. For instance, large numbers of people enter and leave office buildings at approximately the same time. We show that this makes the transformed disturbances an approximately sparse signal, motivating the use of $\ell_1$ regularization.  We test our method via simulations, and results indicate that the method can estimate the thermal dynamic model and transformed disturbance with both open loop and closed loop data, even when the disturbance is not piecewise constant.

	To the best of our knowledge, the only prior work on simultaneously identifying a dynamic model of a building's temperature dynamics and unmeasured disturbances from data are the recent references~\cite{kimcaiaribra:2016,fux2014ekf,hu2016building}. There are many differences between these references and out work. The method proposed in~\cite{kimcaiaribra:2016} estimates the plant parameters and an output disturbance (a disturbance that is added to the plant output) that encapsulates the effect of an unknown input disturbance. In contrast, the proposed method estimates an input disturbance. Both \cite{fux2014ekf} and \cite{hu2016building} take a similar approach: the model is estimated by using data from unoccupied periods (weekends in~\cite{hu2016building}) and assuming that the disturbance is zero during those periods.  Once the model is identified this way, the disturbance is identified using data from occupied periods. Our method uses data collected during regular operation of a building and does not need data collected when the building is empty. Even when data from unoccupied periods is available, assuming the disturbance to be zero during that time is not desirable since doing so will prevent the disturbance from absorbing model mismatch. In contrast to all three methods, the method proposed here can enforce properties of the system that are known from the physics of the thermal processes, such as stability and signs of DC gains for certain input-output pairs. The proposed method consists of solving a convex optimization problem. The methods in \cite{kimcaiaribra:2016,fux2014ekf} require solving non-convex optimization problems. \ifArxivVersion
	The solar heat gain is assumed to be a known input in \cite{Fux:2014}. In
	practice, while solar irradiance can be easily measured with a sensor, there is significant
	uncertainty in the relationship between solar irradiance  and solar heat gain. In the proposed method, we only use solar irradiance as a known input.
	The effective area, which determines the solar heat gain together with solar irradiance, is an not assumed to be known.

	The selection of regularization parameter can be crucial to the identification accuracy of the $\ell_1$ regularized-least-square problems. For the classic Lasso (least absolute shrinkage and selection operator) introduced in \cite{tibshirani1996regression}, there are a number of heuristics. However, the conditions required for ther application do not hold in our case. We therefore propose a distinct heuristic to choose $\lambda$.
	\fi

	The rest of this paper is organized as follows. Section~\ref{sec:problem_formu} formally describes the problem and establishes some properties that will useful later. Section \ref{sec:method} describes the proposed algorithm. We provide evaluation results in Section~\ref{sec:evaluation}. Finally, Section~\ref{sec:conclusion} concludes this work.

	\section{Problem Formulation}\label{sec:problem_formu}
	The indoor zone temperature $T_z$ is affected by three \emph{known}
	inputs: (1) the heat gain added to the zone by the HVAC system,
	\ensuremath{q_\mathrm{hvac}} (kW), (2) the outside air temperature $T_{oa}$ (\degree C), (3) the solar irradiance \ensuremath{\eta^\mathrm{sol}} (kW/m$^2$), and
	the \emph{unknown} disturbance $\ensuremath{q_\mathrm{int}}$ (kW) which is the internal heat gain due to occupants,
	lights, and equipment used by the occupants. So $u(t) \ensuremath{:=}
	[\ensuremath{q_\mathrm{hvac}}(t), T_{oa}(t), \ensuremath{\eta^\mathrm{sol}}(t)]^T \in \mathbb{R}^3$ and $w(t) = \ensuremath{q_\mathrm{int}}(t) \in \mathbb{R}$. The only measurable output is the
	indoor temperature $T_z$, so $y(t) = T_z(t) \in \mathbb{R}$.

	The model we wish to identify is an black box model
	relating the known inputs and the unknown disturbance, to the measured
	output. We will later enforce constraints on the model's parameters by
	relating the model to a physics-based, second-order continuous-time model, making it an ``semi black-box" model.

	\subsection{Discrete-time model to be identified}\label{sec:dt-tf}
	We start with the following 2nd-order discrete-time transfer function model of  the system, with a sampling period $t_s$:
	\begin{align}\label{eq:dtmodel_1}
	y(\ensuremath{z^{-1}}) = &
	\frac{1}{D(\ensuremath{z^{-1}})}\Big[\sum_{j=1}^3[\sum_{i=0}^2\alpha_{ij}z^{-i}]u_j(\ensuremath{z^{-1}})
	\notag \\
	& \quad + [\sum_{i=0}^2\beta_{i}z^{-i}]w(\ensuremath{z^{-1}})\Big]
	\end{align}
	where $D(\ensuremath{z^{-1}}) = 1-\theta_1\ensuremath{z^{-1}}-\theta_2\ensuremath{z^{-2}}$, for some
	parameters $\theta_1,\theta_2$ and $\alpha_{ij},\beta_i$'s, and
	$u[k],w[k],y[k]$ are samples of the continuous-time signals $u(t),w(t),y(t)$. For future
	convenience, we rewrite it as
	\begin{align}\label{eq:dtmodel}
	y(\ensuremath{z^{-1}}) = &
	\frac{1}{D(\ensuremath{z^{-1}})}\Big[K(\ensuremath{z^{-1}})^Tu(\ensuremath{z^{-1}})  + \bar{w}(\ensuremath{z^{-1}}) \Big]
	\end{align}
	where
	\begin{align}\label{eq:K-def}
	K(\ensuremath{z^{-1}}) &:= \begin{bmatrix}
	\theta_3\ensuremath{z^{-2}}+\theta_4\ensuremath{z^{-1}}+\theta_5 \\
	\theta_6\ensuremath{z^{-2}}+\theta_7\ensuremath{z^{-1}}+\theta_8 \\
	\theta_9\ensuremath{z^{-2}}+\theta_{10}\ensuremath{z^{-1}}+\theta_{11}
	\end{bmatrix}
	\end{align}
	and $\bar{w}(\ensuremath{z^{-1}})$ is the z-transform of the \emph{transformed
		disturbance} signal $\bar{w}[k]$ defined as
	\begin{align}\label{eq:def-wbar}
	\bar{w}[k]&:= \beta_0w[k] + \beta_1w[k-1] + \beta_2w[k-2].
	\end{align}
	Performing an inverse z-transformation on
	\eqref{eq:dtmodel}-\eqref{eq:K-def}, yields a difference equation, from which we obtain the linear regression form:
	\begin{equation} \label{eq:regression}
	y[k] = \phi[k]^T\theta, \quad k=3,\dots,\ensuremath{k_\mathrm{max}}
	\end{equation}
	where $\theta^T:=[\theta_p^T,\bar{w}^T]$, in which
	$\theta_p=[\theta_1,\dots,\theta_{11}]^T \in \mathbb{R}^{11}$,
	$\bar{w}=[\bar{w}_3,\dots,\bar{w}_{\ensuremath{k_\mathrm{max}}}]^T \in \mathbb{R}^{\ensuremath{k_\mathrm{max}}-2}$ and
	\begin{multline*
	\phi[k]^T \ensuremath{:=}
	\Big[y[k-1],y[k-2],u_1[k-2],u_1[k-1],u_1[k],\\
	u_2[k-2],\dots,u_2[k],u_3[k-2],\dots,u_3[k], e_{k-2}^T\Big],
	\end{multline*}
	where $e_k $ is the $k$-th canonical basis vector of $\mathbb{R}^{\ensuremath{k_\mathrm{max}}-2}$ in which
	the 1 appears in the $k^\text{th}$ place, and $\ensuremath{k_\mathrm{max}}$ is the number of
	samples. Eq.~\eqref{eq:regression} can be expressed as
	\begin{align}\label{eq:regression-vec}
	y = \Phi\theta,
	\end{align}
	where $y \ensuremath{:=} \left[y[3],\dots,y[\ensuremath{k_\mathrm{max}}]\right]^T \in \mathbb{R}^{\ensuremath{k_\mathrm{max}}-2}$ and
	\begin{align}
	\Phi & \ensuremath{:=}
	\begin{bmatrix}
	\phi[3]^T \\
	\dots\\
	\phi[\ensuremath{k_\mathrm{max}}]^T
	\end{bmatrix} \in \mathbb{R}^{\ensuremath{k_\mathrm{max}}-2 \times \ensuremath{k_\mathrm{max}}+9} \nonumber
	\end{align}
	The problem we seek to address is: \emph{given time traces of inputs and
		outputs, $\{u[k],y[k]\}_1^{\ensuremath{k_\mathrm{max}}}$, determine the unknown parameter vector
		$p \in \mathbb{R}^{11}$ and the unknown transformed disturbance
		vector $\bar{w}:=[\bar{w}_3,\dots,\bar{w}_{\ensuremath{k_\mathrm{max}}}]$, i.e., determine $\theta$}.

	It should be noted that the matrix $\Phi$ is not full column-rank, so there will be an infinite number of solutions to \eqref{eq:regression-vec}. %

	We will use insights from a physics-based model to impose additional constraints on $\theta$.

	\subsection{Insight from an RC network ODE model}\label{sec:RC}
	Figure~\ref{fig:RCschematic} shows a building (left) and a corresponding 2nd-order resistance-capacitance (RC) network model (right). RC-networks are a common modeling paradigm for building thermal dynamics~\cite{Kramer:2012}.

	\begin{figure}[h]
		\includegraphics[width=0.95\linewidth]{2R2C_zone.pdf}
		\centering
		\caption{A photograph for Pugh Hall and a schematic of the ``2R2C'' model.}
		\label{fig:RCschematic}
	\end{figure}
	The ode model of the RC-network model shown in the figure is
	\begin{align} \label{eq:RC}
	\begin{aligned}
	C_z\dot{T}_r &= \frac{T_w-T_z}{R_z}+\ensuremath{q_\mathrm{hvac}}+\ensuremath{A_\mathrm{e}}\ensuremath{\eta^\mathrm{sol}}+\ensuremath{q_\mathrm{int}}\\
	C_w\dot{T}_w &= \frac{T_{oa}-T_w}{R_w}+\frac{T_z-T_w}{R_z}
	\end{aligned}
	\end{align}
	where $C_z, C_w, R_z, R_w$ are the
	thermal capacitances and resistances of the zone and wall,
	respectively, and \ensuremath{A_\mathrm{e}} is the effective area of the building for
	incident solar radiation. All five parameters are positive. Defining the state vector as $x \ensuremath{:=} [T_z,
	T_w]^T \in \mathbb{R}^2$, the ode can be expressed as
	\begin{eqnarray}
	\begin{aligned}\label{eq:ct_ss}
	\dot{x}&=Fx+Gu+Hw, & y &= Jx,
	\end{aligned}
	\end{eqnarray}
	where  $u$,$w$, and $y$ are defined in Section \ref{sec:problem_formu}, and $F\in \mathbb{R}^{2\times 2}, G \in \mathbb{R}^{2 \times 3}$, $H \in \mathbb{R}^{2
		\times 1}$ and $J \in \mathbb{R}^{1 \times 2}$ are appropriate matrices
	that are functions of the parameters $C_z, C_w, R_z,
	R_w,\ensuremath{A_\mathrm{e}}$.
	\ifArxivVersion
	Specifically,
	\begin{align*}
	F & =
	\begin{bmatrix}
	-\frac{1}{C_zR_z} &      \frac{1}{C_zR_z} \\
	\frac{1}{C_wR_z} &     -\frac{1}{C_w}(\frac{1}{R_w}+\frac{1}{R_z})
	\end{bmatrix},\\
	G & =
	\begin{bmatrix}
	\frac{1}{C_z}   &           0             & \frac{A_e}{C_z} \\
	0  &     \frac{1}{C_wR_w}    &     0
	\end{bmatrix}, \\
	H & =
	\begin{bmatrix}
	\frac{1}{C_z}  \\  0
	\end{bmatrix},\\
	J &=
	\begin{bmatrix}
	1 & 0.
	\end{bmatrix}.\\
	\end{align*}
	\fi
	In Laplace domain,
	\begin{align} \label{eq:ct_tf}
	y(s) & = \frac{1}{D(s)} \Big[(s-f_{22}) \left(
	g_{11}u_1(s)+g_{13}u_3(s) \right) \notag \\
	& +f_{12}g_{22}u_2(s) + (s-f_{22})h_{11}w(s) \Big]
	\end{align}
	where $f_{ij},g_{ij},h_{ij}$ are the $i,j$-th entry of the matrices
	$F,G,H$ (respectively) in \eqref{eq:ct_ss}, and
	\begin{align}
	D(s) & = s^2 + d_1s + d_2, \; \text{with} \label{eq:def-D(s)}\\
	d_1 &=
	\frac{1}{C_zR_z}+\frac{1}{C_w}(\frac{1}{R_z}+\frac{1}{R_w}), &
	d_2 &= \frac{1}{C_zC_wR_zR_w}. \notag
	\end{align}
	\emph{We now assume that the discrete-time system~~\eqref{eq:dtmodel_1} was obtained by discretizing the continuous-time system \eqref{eq:ct_tf} by Tustin transformation.} It can be shown through straightforward calculations that the parameters of the discrete-time model -- the $\theta_i$'s -- are related
	to those of the continuous-time model \eqref{eq:ct_tf} as follows:
	\begin{align}
	\begin{aligned}
	&\theta_1 :=\frac{8-2d_2t_s^2}{D_0},
	\theta_2 := -\frac{d_2t_s^2-2d_1t_s+4}{D_0}, \\
	&\begin{bmatrix}
	\theta_3 & \theta_9\\
	\theta_4 & \theta_{10}\\
	\theta_5 & \theta_{11}\\
	\end{bmatrix}  := \frac{t_s}{D_0}
	\begin{bmatrix}
	-2-f_{22}t_s \\-2f_{22}t_s \\ 2-f_{22}t_s
	\end{bmatrix}
	\begin{bmatrix}
	g_{11} & g_{13}
	\end{bmatrix}, \\
	& \begin{bmatrix}
	\theta_6 \\ \theta_7 \\ \theta_8
	\end{bmatrix}  :=
	\begin{bmatrix}
	1 \\ 2 \\ 1
	\end{bmatrix}
	\frac{f_{12}g_{22}t_s^2}{D_0},
	\end{aligned}\label{eq:theta_def}
	\end{align}
	where $D_0=d_2t_s^2+2d_1t_s+4$. Similarly,
	\begin{align}
	[\beta_0,\beta_1,\beta_2] = \frac{t_s\big[(2+\epsilon_0), \;
		2\epsilon_0, \; (-2+\epsilon_0) \big]}{C_zD_0},   \label{eq:def-betas}\\
	\text{ where } \epsilon_0 = -f_{22}t_s =
	\frac{t_s}{C_w}(\frac{1}{R_w}+\frac{1}{R_z}).   \label{eq:def-eps0}
	\end{align}

	\subsubsection{Insight I: Sparsity of transformed  disturbance}\label{sec:sparse_wbar}
	We need a few definitions to talk about \emph{approximately sparse}
	vectors, and \emph{slowly varying} vectors.
	\begin{definition}
		\begin{enumerate}
			\item A vector $x\in \mathbb{R}^{n}$ is $(\epsilon,f)$-sparse if at most
			$f$ fraction of entries of $x$ are \emph{not} in $[-\epsilon,\epsilon]$.
			\item The \emph{change frequency} $c_f(x)$ of a vector $x\in\mathbb{R}^n$ is the fraction
			of entries that are distinct from their previous neighbor: $c_f(x)
			= \frac{1}{n-1}|\{k>1|x_k \neq x_{k-1}\}|$, where $|A|$ denotes the
			cardinality of the set $A$. We say a
			vector $x$ \emph{changes infrequently} if $c_f(x)\ll 1$.
		\end{enumerate}
	\end{definition}
	The following result shows that if the disturbance does not change frequently (which happens if it is piecewise constant), then the transformed disturbance is approximately sparse.
	\begin{proposition}\label{prop:sparse}
		Suppose the disturbance $w[k]$ is uniformly bounded in $k$, it changes infrequently, and $\epsilon_0\ll 1$ where $\epsilon_0$ is defined in \eqref{eq:def-eps0}.  Then, $\bar{w}[k]$ is $(\bar{\epsilon},2c_f(w))$-sparse, where
		$\bar{\epsilon}=\frac{4}{C_zD_0}t_s w_u\epsilon_0$ and $w_u$ is an  upper bound on $|w[k]|$. $\square$
	\end{proposition}

	\begin{proof
		It can be shown from \eqref{eq:def-wbar} and \eqref{eq:def-betas} that
		\begin{align*}
		\bar{w}[k]=&\frac{t_s}{C_zD_0}\big(2(w[k]-w[k-2])\\
		&\quad -\epsilon_0(w[k]+2w[k-1]+w[k-2])\big)
		\end{align*}
		Since $w$ is bounded, $\exists w_l,w_u$ with  $|w_l|\leq w_u$, $w_u\geq0$ s.t. $w[k]\in
		[w_l,w_u]$. Since $c_f(w) \ll 1$ from the hypothesis, for at least  $1-2c_f(w)$ fraction of  $k$'s,
		$w[k]-w[k-2]=0$, and for those $k$'s,
		\begin{eqnarray}
		\begin{aligned}
		\bar{w}[k]&=-\epsilon_0\frac{t_s}{C_zD_0}\big(w[k]+2w[k-1]+w[k-2]\big)\\
		&\in [\frac{-4\epsilon_0t_sw_u}{C_zD_0},\frac{4\epsilon_0t_sw_l}{C_zD_0}]\subset [\frac{-4\epsilon_0t_sw_u}{C_zD_0},\frac{4\epsilon_0t_sw_u}{C_zD_0}],
		\end{aligned}
		\end{eqnarray}
		which proves the result.
	\end{proof}

	Since the product $RC$ is large for large buildings, of the order of few hours~\cite{kimcaiaribra:2016}, it follows from \eqref{eq:def-eps0} that $\epsilon_0$ is small for such buildings. In addition, both $\epsilon_0$ and  $\bar\epsilon$ can be made as small as possible by choosing $t_s$ sufficiently small. The assumption in the proposition,  that $\epsilon_0$ is small, is therefore not a strong one.

	\subsubsection{Insight II:Constraints on parameters}\label{sec:para_const}
	The constraints described below are straightforward to derive, but involve - in a few cases - extensive algebra. We therefore omit the details here; they can be found in the expanded version \cite{ZengSimultaneousArXiV:2017}.
	\paragraph{Stability}
	It can be shown that due to the resistances and capacitances in~\eqref{eq:RC} being positive, the continuous time model~\eqref{eq:ct_tf} is BIBO stable. Since Tustin transformation preserves stability, all poles of the transfer function \eqref{eq:dtmodel_1} should be inside the unit circle~\cite{ogata1995discrete}. It can be shown that this is equivalent to
	\begin{align}
	\begin{aligned}\label{eq:stable_constraint}
	-\theta_2& < 1, & \theta_2+\theta_1&< 1,
	\end{aligned}\\
	\begin{aligned}\label{eq:stable_const_redun}
	\theta_2-\theta_1& < 1.
	\end{aligned}
	\end{align}

	\paragraph{Sign of parameters}
	By using the positivity of the parameters $R_w,R_z,C_w,C_z$, it can be shown
	after some tedious algebra that if $t_s<2min\{\frac{C_wR_wR_z}{R_z+R_w},\sqrt{R_zC_zR_wC_w},\frac{\min(R_zC_z,R_zC_w,R_wC_w)}{3}\}$, the following holds:
	\begin{equation}\label{eq:thetasigns}
	\begin{split}
	&\theta_i>0, i \in \{1,4,5,6,7,8,10,11\}, \\
	&\theta_2<0, \theta_3<0, \theta_9<0
	\end{split}
	\end{equation}
	\paragraph{Positive DC-gain}
	An increase in any of the inputs $\ensuremath{q_\mathrm{hvac}}, T_{oa},
	\eta^{sol}$ represents an increase in the cooling load for the building. A steady state increase in any of these inputs must therefore lead to a steady state increase in the indoor
	temperature $T_z$. In other words, the corresponding DC gains must be positive. Using the previously established fact that the denominator coefficients are positive (see~\eqref{eq:stable_constraint}) it can be shown that positive DC gains are equivalent to
	\begin{eqnarray}
	\begin{aligned}\label{eq:dc_num_const345}
	\theta_3+\theta_4+\theta_5&> 0,
	\end{aligned}\\
	\begin{aligned}\label{eq:dc_num_const678}
	\theta_6+\theta_7+\theta_8&> 0,
	\end{aligned}\\
	\begin{aligned}\label{eq:dc_num_const91011}
	\theta_9+\theta_{10}+\theta_{11}&> 0.
	\end{aligned}
	\end{eqnarray}
	In order to ensure existence of a solution~\cite{luenberger1984linear}, the above constraints are relaxed from a strict inequality to a non-strict one. Additionally, the redundant inequalities~\eqref{eq:stable_const_redun} and~\eqref{eq:dc_num_const678} are removed since they do not change the feasible region,
	\ifArxivVersion
	  where the proof is provided in the Appendix.
	\else
	  where the proof is provided in \cite{ZengSimultaneousArXiV:17}.
	\fi
	The remaining, linearly independent constraints can be compactly written as
	$g(\theta)\leq0$ with $g:\mathbb{R}^{11} \to \mathbb{R}^{15}$, where the inequality is
	entry-wise.

	\section{Proposed method}\label{sec:method}
	Let $S\ensuremath{:=}
	[\bm{0}_{\ensuremath{k_\mathrm{max}}-2\times11} | \;I_{\ensuremath{k_\mathrm{max}}-2}]$ so that
	$S\theta=\bar{w}$. Since we expect $w$ to be piecewise constant and slowly varying, $\bar{w}$ should be approximately sparse (Proposition
	\ref{prop:sparse}). We thus seek a solution to $\mathbf{y}=\Phi\theta$ so that $S\theta$ is sparse, by posing the following optimization problem
	\begin{equation}\label{eq:opt}
	\begin{split}
	\hat{\theta} & = \arg \min_{\theta }   \|y -  \Phi \theta\|_2^2 +\lambda \|S\theta\|_1 \\
	& \quad \text{ s. t. } g(\theta) \leq 0
	\end{split}
	\end{equation}
	where $\lambda > 0$ is a user-defined weighting factor. The $\ell_1$ norm penalty is to encourage sparsity. The estimated plant parameters $\hat{\theta}_p$ and estimated transformed disturbance $\hat{\bar{w}}$ can be recovered from $\hat\theta$ since $\theta^T  = [\theta_p^T, \bar{w}^T]$.

	Regularity of constraints is useful for optimization algorithms to perform well~\cite{Andreani_etal:2014}, and the next result establishes regularity.
	\ifArxivVersion
	The proof is provided in the Appendix.
	\else
	The proof of the result is straightforward and available in the expanded version \cite{ZengSimultaneousArXiV:17}, but omitted here due to space constraints.
	\fi

	We call a point $\theta$ \emph{physically meaningful} if none of the
	three SISO transfer functions in~\eqref{eq:dtmodel} is identically zero.
	\begin{proposition}\label{prop:feasible_convex_regular}
		The optimization problem \eqref{eq:opt} is feasible, convex, and every physically meaningful feasible $\theta$ is a regular point of the constraints.
	\end{proposition}

    The optimization problem \eqref{eq:opt} is convex. All numerical results presented in this paper were obtained by using the \texttt{cvx} package for solving convex problems in MATLAB$^\copyright$ ~\cite{cvx}.

	\subsection{Regularization Parameter Selection}\label{sec:regu_select}
	\ifArxivVersion %
	In contrast to the classic Lasso problem, which is to minimize $\|\bm{y} - \Phi\theta \|_2^2$+$\lambda\|\theta\|_1$~\cite{tibshirani1996regression}, our problem falls into the category of ``generalized Lasso'' due to the replacement of $\theta$ by $S\theta$ in the cost function~\cite{duan2016generalized}. Two common heuristics for choosing $\lambda$ for the classic Lasso problem are cross-validation~\cite{pendse2011tutorial} and L-curve-based curvature methods~\cite{hansen1992analysis}. However, neither of them is applicable to our problem as they have implementation requirements that our problem does not satisfy.

	Cross validation divides datasets into $K$ folders and requires that parameters to be retrieved are the same among such folders, whereas parameters in our problem contain transformed disturbance, which may differ from one day to the next. The L-curve, which is ``a log-log plot of the norm of a regularized solution $\|\theta\|_1$ versus the norm of the corresponding residual norm $\|y - \Phi \theta\|_2$'', can graphically display the trade-off between the size of a regularized solution and its fit to the given data, where optimal regularization parameter that minimize the trade-off lies at the corner of such L-curve. For the $L$-curve method, a solution path that changes monotonically with respect to $\lambda$ is essential, i.e., $(\Phi^T\Phi)^{-1}$ needs to be diagonally dominant~\cite{duan2016generalized}. That too is also not satisfied in our case.
	\else
	The problem \eqref{eq:opt} can be converted to the classic Lasso, which is to minimize $\|\bm{y} - \Phi\theta \|_2^2$+$\lambda\|\theta\|_1$~\cite{tibshirani1996regression}. There are a number of heuristics for choosing $\lambda$ in the Lasso, such as L-curve \cite{hansen1992analysis} and cross-validation~\cite{pendse2011tutorial}. However, the conditions for application of these heuristics are not satisfied in our case \cite{ZengSimultaneousArXiV:17}.
	\fi

	We use the following heuristic to choose $\lambda$, which is inspired by the L-curve method. First, plot both the solution norm and residual norm individually against $\lambda$ by repeatedly solving Problem~\eqref{eq:opt}  for various values of $\lambda$.
	An illustration of these two plots is shown in Figure \ref{fig:L1-L2-curve}.
	Second, identify a value $\lambda_1$ so that the solution norm is smaller than a user-defined threshold for any $\lambda>\lambda_1$, and then identify $\lambda_2$ so that the residual norm is smaller than a user-defined threshold for any $\lambda<\lambda_2$. If $\lambda_2>\lambda_1$, choose $\lambda$ to be $\lambda_1$. If not, pick another threshold, and continue until this condition is met.

	\begin{figure}[htb]
		\includegraphics[width=0.8\linewidth]{L1_L2.pdf}
		\centering
		\caption{Illustration of regularization parameter selection}
		\label{fig:L1-L2-curve}
		\vspace{-0.2 cm}
	\end{figure}

	\section{Evaluation of Proposed Algorithm}\label{sec:evaluation}
	The continuous-time RC model~\eqref{eq:ct_ss}  is used to generate training and
	validation data. The parametrs of the model were chosen by manual calibration of the model to data collected from a large auditorium in a campus building at the
	University of Florida (Pugh Hall; shown in Figure \ref{fig:RCschematic}). Four scenarios are tested:
	\begin{enumerate}
		\item \emph{OL-PW:} Open-loop with piecewise-constant disturbance;
		\item \emph{OL-NPW:} Open-loop with not piecewise-constant disturbance;
		\item \emph{CL-PW} Closed-loop with piecewise-constant disturbance;
		\item \emph{CL-NPW:} Closed-loop with not piecewise-constant disturbance;
	\end{enumerate}
	If the disturbance  $w[k]$ is piecewise constant, since that is slowly varying, the transformed disturbance $\bar{w}[k]$ will be approximately sparse.
	\ifArxivVersion
	The algorithm is expected to perform well in the OL-PW scenario since it satisfies the piecewise constant assumption the method is based on, and identification with open-loop data is generally easier than with closed loop~\cite{ljungbook:1999}. The CL-NPW scenario is the most relevant in practice, but it is likely to be the most challenging for the method.
	\fi
	%
	In all four scenarios, the same
	input data sequence for ambient temperature (from \url{weatherunderground.com})
	and solar irradiance data (from NSRDB: \url{https://nsrdb.nrel.gov/}), both for Gainesville, FL, are used. In the two open-loop scenarios, the input component $\ensuremath{q_\mathrm{hvac}}$ is	somewhat arbitraruly chosen, while in the two closed-loop scenarios,
	it is decided by a PI-controller that tries to maintain the zone temperature near a setpoint $T^\mathrm{ref}$. To have exciting input to aid in identification, the setpoint $T^\mathrm{ref}$ is chosen to be a PRBS
	sequence~\cite{ljungbook:1999}. To ensure that
	occupant comfort is not compromised, the setpoint is constrained to
	lie within $22$\degree C and $27$\degree C.
	The disturbance signal $\ensuremath{q_\mathrm{int}}$ is picked somewhat arbitrarily during manual calibration of the RC network model to Pugh Hall data. The training data are shown in Figure~\ref{fig:week_input}.
	\begin{figure}[htb]
		\centering
		\includegraphics[width=0.95\linewidth]{inputs_output.pdf}
		\caption{Training data for algorithm evaluation. The data  $\ensuremath{\eta^\mathrm{sol}}, T_{oa},\ensuremath{q_\mathrm{int}}$ shown here are used in all four scenarios;
			$\ensuremath{q_\mathrm{hvac}}, T_z$ shown here are for the CL-NPW scenario.}
		\label{fig:week_input}
	\end{figure}
	Notice from the figure that the disturbance $\ensuremath{q_\mathrm{int}}$ is large; sometimes as large as the cooling power provided by the HVAC system.

	\subsection{Plant identification results}\label{Sec:results-plantid}
	\paragraph{Parameters}
	Table~\ref{tab:theta-compare} shows the true values of the plant
	parameters, $\theta_1,\dots,\theta_{11}$. It also shows the
	corresponding estimation errors (in percentage) for the OL-PW and
	CL-NPW scenarios. We can see from the table that performance of
	the method is similar with both open-loop and closed-loop data.  Second, the two parameters that determine the characteristic equation
	are estimated highly accurately. Third, there is more error in the  estimate of numerators. While some are more accurate than others,
	the numerator coefficients corresponding to the input $\ensuremath{\eta^\mathrm{sol}}$ has the most error. A possible reason for this high error is the lack of richness in the $\ensuremath{\eta^\mathrm{sol}}$ data. See Figure \ref{fig:week_input}: $\ensuremath{\eta^\mathrm{sol}}$ is the least rich among all the input signals.

	Results for the remaining two scenarios are similar, but are not shown due to space constraints.
	\begin{table}[t]
		\centering
		\caption{Plant parameters and errors in their estimates.}
		\label{tab:theta-compare}
		\begin{tabular}{c|c|c|c|c}
			\multicolumn{2}{c|}{\multirow{2}{*}{$\theta$}} & \multicolumn{2}{c|}{$\frac{\theta-\hat{\theta}}{\theta}\%$} &
			\multirow{2}{*}{input}\\
			\multicolumn{2}{c|}{}                  &  (OL-PW)   &  (CL-NPW)   &  \\ \hline
			$\theta_1   $ & $ 1.98\times 10^{-0} $ & $  -0.075$ & $   0.042 $ &        \\
			$\theta_2   $ & $-9.76\times 10^{-1} $ & $  -0.151$ & $   0.085 $ &        \\ \hline
			$\theta_3   $ & $-4.35\times 10^{-3} $ & $  -9.214$ & $  -8.024 $ &        \\
			$\theta_4   $ & $ 5.21\times 10^{-5} $ & $ -59.48 $ & $-108.2   $ & \ensuremath{q_\mathrm{hvac}} \\
			$\theta_5   $ & $ 4.40\times 10^{-3} $ & $  -7.493$ & $  -6.36  $ &        \\ \hline
			$\theta_6   $ & $ 1.86\times 10^{-5} $ & $ -18.64 $ & $ -48.90  $ &        \\
			$\theta_7   $ & $ 3.72\times 10^{-5} $ & $  38.15 $ & $  22.35  $ & $T_{oa}$  \\
			$\theta_8   $ & $ 1.86\times 10^{-5} $ & $ -39.89 $ & $ -68.32  $ &        \\\hline
			$\theta_9   $ & $-3.05\times 10^{-2} $ & $-112.6  $ & $-232.1   $ &        \\
			$\theta_{10}$ & $ 3.65\times 10^{-4} $ & $-12300  $ & $-19320   $ &\ensuremath{\eta^\mathrm{sol}} \\
			$\theta_{11}$ & $ 3.08\times 10^{-2} $ & $  33.18 $ &   $-2.881 $ &
			\\ \hline
		\end{tabular}
	\end{table}

	\paragraph{Frequency response} For prediction accuracy, frequency response is more important than individual parameters. %
	\ifArxivVersion
	Figure~\ref{fig:bode_OL_cmp} compares the frequency response of the
	identified plants with their true values for the two open loop
	scenarios. Notice that just as in case of parameters, the estimates
	corresponding to the input $\ensuremath{\eta^\mathrm{sol}}$ is the poorest. We believe this
	is due to the lack of sufficient excitation in the data;
	cf. Figure~\ref{fig:week_input}. The input $T_{oa}$ also has low excitation at higher
	frequencies, and therefore has poor estimates in higher frequencies.

	For the transfer functions from inputs $\ensuremath{q_\mathrm{hvac}}$ to output $T_z$, the maximum absolute error in the estimated frequency response is:
	\begin{align*}
	\max_{\omega}\frac{|\hat{G}_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega) - G_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega)|}{|G_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega) |} = 0.108
	\end{align*}
	and occurs at $ \omega = 1/(10$ weeks). The maximum errors for the transfers functions from $T_{oa}$ and $\eta_{sol}$ to $T_z$ occur at the Nyquist frequency.
	\begin{figure*}[htb]
		\includegraphics[width=0.95\linewidth]{bode_OL_cmp.pdf}
		\centering
		\caption{Bode magnitude plots of true and identified systems for OL scenarios.}
		\label{fig:bode_OL_cmp}
	\end{figure*}
	\begin{figure*}[htb]
		\includegraphics[width=0.95\linewidth]{bode_CL_cmp.pdf}
		\centering
		\caption{Bode magnitude plots of true and identified systems for CL scenarios.}
		\label{fig:bode_CL_cmp}
	\end{figure*}
	\else

	Fig.~\ref{fig:bode_jnl_cmp} shows the Bode plots of
	the true and identified models. Due to space limitations, we only show the Bode plots for the OL-PW and CL-NPW  scenarios.
	\begin{figure}[htb]
		\includegraphics[width=1\linewidth]{figures/bode_jnl_cmp.pdf}
		\centering
		\caption{Bode magnitude plots of the true and identified systems.}
		\label{fig:bode_jnl_cmp}
	\end{figure}
	For the transfer functions from inputs $\ensuremath{q_\mathrm{hvac}}$ to output $T_z$, the maximum absolute error in the estimated frequency response is:
	\begin{align*}
	\max_{\omega}\frac{|\hat{G}_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega) - G_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega)|}{|G_{\ensuremath{q_\mathrm{hvac}} T_z}(j\omega) |} = 0.108
	\end{align*}
	and occurs at $ \omega = 1/(10$ weeks). The maximum errors for the transfers functions from $T_{oa}$ and $\eta_{sol}$ to $T_z$ occur at the Nyquist frequency.

	\fi

	\subsection{Disturbance}
	The estimated transformed disturbance, $\hat{\bar w}$, for all four scenarios are shown in Fig.~\ref{fig:Wbar_cmp}. The estimates are quite accurate when the true values are large, but less accurate otherwise. However, the estimates capture the trend of the true values quite accurately, even when the true disturbance is not piecewise constant, in which case the transformed disturbance may be neither approximately sparse nor slowly varying.

	\ifArxivVersion
	\begin{figure*}[htb]
		\includegraphics[width=1\linewidth]{wbar_results_cmp.pdf}
		\centering
		\caption{Comparison of identified and actual transformed disturbance. Bottom two plots are zoomed version on Tuesday of the top two plots.}
		\label{fig:Wbar_cmp}
	\end{figure*}
	\else
	\begin{figure}[htb]
		\includegraphics[width=1\linewidth]{figures/wbar_results_cmp.pdf}
		\centering
		\caption{Comparison of identified and actual transformed disturbance. Bottom two plots are zoomed version on Tuesday of the top two plots.}
		\label{fig:Wbar_cmp}
	\end{figure}
	\fi
	\subsection{Validation through temperature prediction}
	The plant identified with data from one week is used to predict temperatures in another week. The disturbance data is the same between the calibration and validation data sets but the input $u$ and output
	$y$ data sets are distinct. The rms value of the prediction error of zone temperature is 1.2 \degree C for OL-PW case and 0.1 \degree C for CL-NPW case; see Figure \ref{fig:Tz_cmp}.
	\begin{figure}[htb]
		\includegraphics[width=0.9\linewidth]{Tr_results_cmp.pdf}
		\centering
		\caption{Predicted and actual zone temperature (validation dataset). Upper: OL-PW; Lower: CL-NPW.}
		\label{fig:Tz_cmp}
	\end{figure}
	\ifArxivVersion
	For validation, we use the solar irradiance and ambient temperature from the validation data set; the temperature setpoint is arbitrarily picked as another PRBS sequence lies within  $[22,27]$ \degree C, and the disturbance is piecewise-constant (see Figure~\ref{fig:week_input}). Compared to the large inaccuracies in the estimated plant parameters, predictions of the zone temperature, shown in Figure~\ref{fig:Tz_cmp}, are much smaller. As we can see from the figure, the error is more pronounced in some days of the week, while extremely small in other days. \fi %

	\section{Conclusion}\label{sec:conclusion}
	The main advantages of the method is posing the estimation problem as a convex optimization problem with constraints from physical insights about the system and the disturbances, without requiring specially collected data. Previous methods lacked both convexity and/or physically meaningful constraints. The main limitation is that the identified disturbance is a linear transformation of the true disturbance with unknown coefficients. This presents a challenge in verifying the disturbance estimates when the method is applied to data from a real building. Extracting $w$ from $\bar{w}$, so that the estimate can be verified in a test setting, is a topic of future work.
	\bibliographystyle{IEEEtran}