A Transfer Operator Methodology for Optimal Sensor Placement Accounting for Uncertainty
Himanshu Sharma Umesh Vaidya Baskar Ganapathysubramanian









mymainaddress]Himanshu Sharma


[mycorrespondingauthor]Corresponding author
mysecondaryaddress]Umesh Vaidya


mymainaddress]Baskar Ganapathysubramanianmycorrespondingauthor

baskarg@iastate.edu

[mymainaddress]Department of Mechanical Engineering, Iowa State University, Ames, IA 50010, USA

[mysecondaryaddress]Department of Electrical and Computer Engineering, Iowa State University,Ames,IA,50010,USA










Sensors in buildings are used for a wide variety of applications such as monitoring air quality, contaminants, indoor temperature, and relative humidity. These are used for accessing and ensuring indoor air quality, and also for ensuring safety in the event of chemical and biological attacks. It follows that optimal placement of sensors become important to accurately monitor contaminant levels in the indoor environment. 
However, contaminant transport inside the indoor environment is governed by the indoor flow conditions which are affected by various uncertainties associated with the building systems including occupancy and boundary fluxes. Therefore, it is important to account for all associated uncertainties while designing the sensor layout. 
The transfer operator based framework provides an effective way to identify optimal placement of sensors. Previous work has been limited to sensor placements under deterministic scenarios. In this work we extend the transfer operator based approach for optimal sensor placement while accounting for building systems uncertainties. The methodology provides a probabilistic metric to gauge coverage under uncertain conditions. We illustrate the capabilities of the framework with examples exhibiting boundary flux uncertainty. 



CFD, Uncertainty, Optimal Sensor placements, Perron-Forbenious Operator.














§ INTRODUCTION

With increasing miniaturization in electronic components the deployment of multiple sensors is becoming an integral part of daily life. These sensors are typically interconnected (either locally or centrally), and make up a sensor network. The spatial placement of sensors that make up the network is especially important in the context of ensuring indoor air quality (IAQ) in the built environment. This is because an individual in the western world spends about 90% of their lifetime indoor (World Health Organization <cit.>). Degraded IAQ is linked to various respiratory diseases <cit.>, with long term exposure linked to serious illness and death <cit.>. Additionally, in large public spaces such as atriums, airports, subway stations, and indoor stadiums the assessment of indoor air quality is critical from the security and safety standpoint. This is especially important to avoid transmission of infectious diseases (TID) such as influenza, tuberculosis, and SARS. The spread of such diseases in the built environment has been widely studied <cit.> and is closely connected with air flow and contaminant/infection  transport. An optimally placed set of sensors can effectively mitigate these problems. 

For designing sensor networks to monitor IAQ in indoor spaces, a complete knowledge of air-flow distribution is needed. Additionally, the sensor network must satisfy several criteria such as coverage, sensitivity, response time and detection threshold <cit.>. Broadly, there are three approaches for sensor placement in the built environment: 1) engineering and heuristic methods, 2) optimization and inverse methods, and 3) dynamical systems based methods. The engineering methods are mostly thumb rule based. They usually do not take into account the flow in the space. This affects the total coverage of the sensor network and leaves some regions uncovered. While very simple to use, this approach results in expensive sensor networks, as more number of sensor are required to improve the coverage of the space <cit.>. The optimization and inverse methods were developed to overcome the shortcomings of the heuristic approach. These methods use the actual flow fields or flow rates – either from CFD simulations or from zonal models. A cost function is constructed that encodes the coverage of sensing along with additional constraints such as cost, and detection threshold. This cost function is then minimized to identify the optimal sensor locations <cit.>.

The optimization based techniques has been successfully used for a wide variety of sensor placement scenarios <cit.> using both CFD or zonal simulations. While these methods account for the effects of flow fields in the sensor placement process, challenges include computational difficulties in finding the optima. Furthermore, these methods rely on solving the advection-diffusion equation for multiple release scenarios and are based on iterative solves which makes them very computationally challenging. Similar issues exist with the inverse and adjoint based methods <cit.>. The choice of the initial guess is extremely crucial in these approaches, and typically work well when the number of sensors to be placed is small. Finally, this approach can only find locally optimal sensor locations, with no guarantees of global optima.

The third category of methods is based on a dynamical systems approach to identify the sensor placement locations <cit.>. The approach is based on the construction of the Perron-Frobenius (PF) operator to determine the advection of a contaminant subjected to a non-linear flow field. The key advantage behind the PF operator framework is that the nonlinear evolution of contaminant under nonlinear and uncertain flow fields can be replaced by a equivalent (but high-dimensional) linear evolution of contaminant density. The linear nature of the PF operator affords three significant capabilities: First, this allows well developed concepts from linear control theory to be applied to contaminant transport analysis, second, the resulting algorithms are based on simple matrix-vector products, which make implementation and deployment extremely simple and efficient (both in terms of memory as well as computation), and third, the approach provides guarantees of finding the globally optimal sensor locations. Fontanini et. al <cit.> showed the construction and usage of the discrete PF operator for fast computation of contaminant transport in the built environment. They then used the PF operator to identify optimal sensor locations under deterministic flow scenarios <cit.>. 

In recent years there has been increasing interest in incorporating the effects of various uncertainties in the building during analysis of building performance <cit.>. The uncertainties associated with buildings can be broadly classified into three categories: (1) environmental, (2) workmanship and quality of building elements, and (3) behavioural. Environmental factors consist of location-specific weather conditions which a building is subjected to. This appears as uncertain boundary conditions. Workmanship and quality is associated with the quality and design of building and construction materials. This appears as uncertainty in material properties used in the analysis. Finally, the uncertainties associated with human behaviour, such as occupancy, interior design and usage of appliances are categorized under behavioural uncertainty. Previous studies have shown the importance of these uncertainty on the energy demands <cit.> as well as on the indoor air quality <cit.>. Such prior work has emphasized the importance of accounting for the uncertainty associated with building systems, and motivates this study to develop a sensor placement strategy which can account for uncertainties. 

In the present work we extend the dynamical systems based approach developed by Fontanini et.al <cit.> to account for various building uncertainties while performing sensor placements. We specifically focus on uncertainties associated with environmental effects as they have been shown to significantly affect the flow dynamics within the built environment[The occupancy affects the flow dynamics, and we have previously shown how to account for this uncertainty <cit.>. The quality and workmanship does plays a relatively smaller role due to the fact that there are well defined standards (like ASTM,ASHRAE) that minimize variability]. Our main contributions are as follows: (a) we contribute to the development of robust sensor placement methods by developing a novel optimal sensor placement algorithm which can account for environmental  uncertainties and can be easily extended to account for other uncertainties like occupant behavioural uncertainties, (b) we show implementation of the approach for actual complex building geometries, and (c) show incorporation of sensing constraints such as sensor location, sensing area and accounting for sensor quality to plan the sensor network.

The outline of the paper is as follows: We first briefly outline the deterministic approach in Section <ref> before extending the formulation to include uncertain effects in Section <ref>. In Section <ref> we illustrate the method with representative examples using 2D and 3D complex geometries. Finally, we detail conclusions and future avenues of work in Section <ref>.



§ BRIEF OVERVIEW OF THE DETERMINISTIC PF APPROACH



 §.§ Problem definition

Consider a domain Ω⊂ℝ^d and a given non-linear (un)steady flow field, 𝐔, in the domain Ω. Suppose a contaminant is released in the domain, the aim is to find the sensor locations in Ω that will optimally detect the contaminant release, within the response time of sensors. The sensor placement problem can be posed in a few different ways, depending on what constraints are important to respect. Different constraints result in the following scenarios for sensor placement: (a) Placement when the number of sensors is fixed; (b) Placement when a threshold coverage is desired, with no limit on the number of sensors; and (c) Placement when there are sensing and placement constraints. The third scenario can be further divided into cases where (i) the sensor locations are constrained to be placed outside the occupied zone, (ii) the sensor coverage is only needed in the occupied zone, and (iii) sensor locations are constrained to be placed outside the occupied space, and sensor coverage is needed only for the occupied space. We briefly review the PF based approach developed in Fontanini et.al <cit.> that naturally accounts for all these scenarios.



 §.§ A brief overview of the deterministic approach for PF based sensor placement
  
The dynamical systems approach  <cit.> utilizes the Perron-Frobenious (PF) operator for modeling the evolution of contaminants under non-linear vector fields. The contaminant history is then used to formulated the problem of sensor placement by accounting for various constraints. We first briefly describe how the PF operator is constructed, followed by the deterministic sensor placement approach. Following this, we extend to the uncertain case.  



  §.§.§ Contaminant Transport and Construction of the PF Operator




Contaminant transport is modelled by the advection-diffusion partial differential equation given by Eq-<ref>. 
    ∂Φ/∂ t + ∇(UΦ) +∇^2(DΦ) = S_Φ

The contaminant density, Φ, is propagated by the air flow field, U. The flow field can be generated either experimentally, or computationally using computational fluid dynamics (CFD). Here D is the diffusion constant and S_Φ is the source term. The velocity flow field, U, can be steady or can be a function of time. Eq.<ref> can be thought of as an operator (acting on the field Φ) transporting contaminant from time t to t+δ t  <cit.>. This operator is called the Perron-Frobenius (PF) operator L(·) Eq.(<ref>). 

    Φ_t+δ t = L(Φ_t)

The discrete form of this operator is called the Markov matrix 𝐏. The Markov matrix is a non-negative square matrix with several useful properties including row-, and column- stochasticity (i.e. row- and column- elements sum to 1) and unit eigen spectrum. There are several approaches to calculate the entries of the Markov matrix, especially for contaminant transport <cit.>. We specifically use the Eulerian based method <cit.>, as this places no restrictions on the Markov time step <cit.>.

Once constructed, the Markov matrix can be used to propagate the concentration field at discrete time instances, t_i, to reach a time horizon t_f = t + m δ t. This propagation is performed by simple matrix-vector multiplication [We abuse notation to interchangeably refer to Φ as a field, and as its equivalent vector form], Eq.(<ref>). 

    Φ_t_i+1 = Φ_t_i𝐏 + Ŝ_t_i,t_i+1 i ∈{0,…,m }

The source term, Ŝ_t_i,t_i+1, includes volumetric and inlet sources in the domain <cit.>. An example demonstrating the use of this operator is shown in Fig.<ref>. 

The next step after constructing the 𝐏 matrix is to calculate the contaminant tracking matrix 𝐐_τ.




  §.§.§ Contaminant tracking matrix

 The Markov matrix entries correspond to the probabilities of moving from one state to another in a time step δ t. The next step is to use this matrix to construct the contaminant tracking matrix 𝐐_τ, to find how a constant contaminant source is transported over time interval τ = m×δ t. This can can be done as follows:

    𝐐_τ =  𝐈 + 𝐏 + 𝐏^2 + 𝐏^3 + … + 𝐏^m

𝐐_τ is the matrix which encodes the contaminant tracking history, i.e. where the contaminant will propagate in time τ. Figure-<ref> shows how we construct the contaminant transport for a representative problem of airflow in an aircraft cabin studied by Chen et. al <cit.>.


Once the contaminant tracking matrix is constructed the sensor locations can be decided using 𝐐_τ. 



  §.§.§ Sensor placement

The contaminant tracking matrix, 𝐐_τ, encodes all release scenarios. Each row, i, of the matrix encodes a release scenario, with contaminant release from state i. The non-zero elements of that row indicate which states contain contaminant that was released from state i. It was shown by Vaidya et. al. <cit.> that the optimal states to place sensors are those that are most observable, i.e. the columns of 𝐐_τ that are most non-zero. Choosing the sensor locations becomes the so-called set cover problem of combinatorics <cit.>. This essentially identifies a set of locations that maximize the coverage of states. A greedy algorithm <cit.> is used to find the sensor locations according to the maximum number of covered states.




  §.§.§ Applying sensing constraints
 
There are various constraints associated with sensors that have to be accounted for while deciding on the sensor location. Each sensor has an associated measurement accuracy threshold, which depends on the type/quality of the sensor. The accuracy threshold can be accounted for by inspecting the column entries of the contaminant transport matrix, 𝐐_τ. The entries with larger values represent a stronger signal. The objective of this step is to replace those entries from the contaminant transport matrix which are below a sensing threshold. The threshold value is a non-dimensional sensing threshold that is the ratio of value detected to the product of the source release rate S_source and the release time τ, and is defined as
ϵ_acc = μ_detect/S_sourceτ = μ_detect/μ_source. The values ϵ_acc are usually prescribed by the sensor manufacturer. Since we have constructed the matrix 𝐐_τ we can apply this threshold on the matrix as Eq.(<ref>).

    𝐐^*_τ = 
     1    if 𝐐_τ≥ϵ_acc
    
    0,                 otherwise

This equation serves as an operator that converts a matrix with real entries, 𝐐_τ, to a matrix with binary entries, 𝐐^*_τ. The entries in 𝐐^*_τ that are 1 correspond to the states that can be sensed by sensors with accuracy of ϵ_acc if a sensor is located in the state corresponding to the column in which the entry resides.

Analogously, not every location in the domain is suitable to place a sensor due to practical reasons such as aesthetics, sensor installation limitations, and operation requirements (occupants can affect the sensor operation by stepping over them). Therefore, it is important to account for constraints associated with sensor placement selection. In a broad sense these constraints can classified as:  1) location constraint, where the sensors must be placed outside the occupied zone of the building; and 2) sensing constraint, where we are interested in only the occupied zone. To account for the location constraint, if (ℕ^nloc) states in the domain are not able to accommodate a sensor then the columns corresponding to those states are removed and replaced by zeros: 

    𝐐^*_τ(:,j) = 0 ;      j ∈ℕ^nloc

Similarly, for accounting for the sensing constraints, the states (ℕ^nsen) which are not in the occupied space are unimportant, and are replaced by zeros:

    𝐐^*_τ(i,:) = 0 ;      i ∈ℕ^nsen






§ EXTENSION TO THE UNCERTAIN CASE

We have detailed the PF based methodology for sensor placement under a deterministic flow field. A building in general can be subjected to various weather scenarios leading to different HVAC conditions, interior arrangements and occupancy conditions. These changes brings uncertainty in the velocity field (airflow) inside the building. Therefore, it is important to account this distribution of velocity fields that can occur to identify the sensor locations which can perform optimally in every flow condition. Wang et.al <cit.> showed that HVAC performance is significantly affected by weather conditions. The HVAC operation and weather conditions together significantly affect the flow dynamics in the building. Therefore, in the current study we illustrate the sensor placement approach by accounting for weather uncertainty. 



 §.§ Problem definition under uncertainty


Consider the domain Ω⊂ℝ^d with the boundary ∂Ω. 
We consider a set ξ(t) ∈ℝ^n of random variables, where n is the number of random variables which affect the flow field [We assume that all random variables are independent identically distributed (i.i.d.). Note that at places we ignore the temporal dependence of the random variable ξ for notational convenience.]

As discussed earlier, the flow field is affected by the boundary conditions (environment) of the domain, occupancy and material properties. These uncertainties can be represented in terms of a set of random variables, ξ <cit.>. 


The flow field in the domain is now a function of these random variables, 𝐔≡𝐔(ξ(t)). This definition naturally allows extension to account for other uncertain parameters.  With 𝐔(ξ(t)), stochastic, the advection-diffusion equation (<ref>) transforms to a stochastic partial differential equation. 

The problem definition now becomes: find sensor locations that optimally sense contaminant distributions under uncertain flow conditions.  In contrast to the deterministic case that maximizes the coverage using the PF operator for the deterministic flow field, we will use PF operator constructed for the stochastic flow field to maximize the coverage. We show that this approach essentially lead to the maximization of the expectation of the coverage.

The construction of PF operator for the case when the velocity field, 𝐔(ξ), is stochastic will proceed as follows: We sample the stochastic velocity field and create a finite set of possible realizations with associated probabilities of occurrence. We then construct the PF operators for each of these realizations, and construct the contaminant transport matrix, 𝐐_τ. We then evaluate the expected coverage using the set of contaminant transport matrices and their associated probabilities. We detail each of these steps next. 



 §.§ Constructing PF operators and Contaminant Transport for the sampling points


We construct a finite set of realizations of a random variable by choosing M samples from the distribution of ξ. The M samples forming the set 𝐒 = {s_1,…,s_M} are chosen such that they represent the ξ distribution in a statistical sense. Specifically, associated with each sample is a probability of occurrence (weight) represented as Θ={θ_1,…,θ_M} satisfying, θ_i≥ 0 and  ∑_i θ_i=1. Thus, each sample, s_i has an associated flow field, 𝐔(s_i), along with a probability of occurrence of this flow field. 

The set, (𝐔(𝐒), Θ) is chosen such that they produce the same expectation as the full distribution,

    ∫𝐔(ξ) dξ = ∑_i=1^Mθ_i*𝐔(s_i)

Given the probability distribution of the random variable there are standard approaches to choosing the sampling points and associated weights <cit.>. These are discussed in detail in appendix-<ref>.


In the current work we use CFD simulations to calculate the flow field corresponding to each sample point. In our case, the sample points correspond to realizations of different boundary conditions arising from weather variability that affect the flow field. Therefore, for each sample point a flow field is computed which is collectively represented as 𝐔 = {𝐔_1,…,𝐔_M}. Now, using the approach discussed in section-<ref> we construct the PF operators 𝐏 = {𝐏_1,…,𝐏_M} corresponding to each velocity field in 𝐔. Under the i.i.d. assumption on the random variable ξ(t), the stochastic counterpart of the deterministic PF operator (Eq. <ref>) for the propagation of contaminant can be written as <cit.>

    Φ_t_j+1 = Φ_t_j∑_i=1^M θ_i 𝐏_i + Ŝ_t_j,t_j+1 j ∈{0,…,m }

This equation indicates that the PF operator for the stochastic flow field is just the probability weighted sum of the individual PF operators i.e., the expected PF operator. Now due to the linear nature of the proposed transfer operator framework, the use of expected PF operator,  ∑_i=1^M θ_i 𝐏_i, to maximize the coverage will essentially lead to maximization of the expectation of the coverage.

Next, the contaminant transport matrices 𝐐={𝐐_τ,1,…,𝐐_τ,M} for each element of the set 𝐔 are constructed. The final step is to threshold the set of matrices based on sensor accuracy, and apply any location and sensing constraints. The procedure discussed in section-<ref> can be applied to each contaminant tracking matrix by applying Eqs. <ref>,<ref>,<ref>, resulting in 𝐐^*_τ={𝐐^*_τ,1,…,𝐐^*_τ,M}. 





 §.§ Calculation Expected Volumetric Coverage and Sensor locations

Coverage for a sensor network is interpreted as how well the sensors can monitor the domain. This is quantified in terms of what fraction of the volume of the domain is being sensed. Considering the fact that the CFD simulations are usually performed on a nonuniform discretization of the domain, we scale the value of each state in 𝐐^*_τ,i with the volume of that state. Thus 𝐐^** is the volumetric scaled version of 𝐐^*, and each row of it is computed as

    𝐐^**(i,:) = 𝐐^*(i,:)*V_ω_i/ V_Ω tot

where V_ω_i is the volume of each state and V_Ω tot total volume of domain Ω. We compute the volumetrically scaled version of the coverage for each contaminant tracking matrix to get 𝐐^** =  {𝐐^**_1,⋯, 𝐐^**_M}. Next, we calculate the total coverage produced if a sensor is placed in state j, for every state j = 1,⋯, N in the domain. This is constructed simply by computing the column sum of the 𝐐^** matrices

    𝐯̅(j) = ∑_k=1^N 𝐐^**(j,k)

This operation is mathematically equivalent to the column wise L_1 norm of the matrix, and results in a vector 𝐯̅ of size N × 1. We perform this operation to get the coverage for each of the M samples to get

    𝕍 = {𝐯̅_1, 
        …,𝐯̅_M}




Finally, the expected coverage is evaluated as the weighted sum of the coverage of each sample, 𝐯̅_M, with the associated probabilities θ_M. The expected coverage, 𝔼[𝕍], is obtained as

    𝔼[𝕍] = ∑_i=1^M𝐯̅_̅i̅θ_i

Once we have calculated the expected coverage vector, 𝔼[𝕍], the location of the first sensor,  k̅(1), is the state (i.e. element of the expected coverage vector), k of 𝔼[𝕍] with maximum value:

    k̅(1) = _k 𝔼[𝕍̅](k)

As in the deterministic case, the identification of the next sensor is done iteratively after removing the coverage of the previously placed sensors <cit.>.  
To find the next sensor, the information of the first sensor is removed from consideration (since we want to ensure that the next sensor has coverage over parts of the domain that are not covered by the first sensor). This is done by zeroing out the row and column corresponding to the state, k̅(1), for each  (𝐐^**)_i matrix. Then the sensor placement procedure is repeated ( i.e. Eq.<ref>-<ref>) until all the sensors are placed in the domain or the required criteria for coverage is full filled. 



 §.§ Overview of algorithm

 The complete expectation based algorithm is illustrated as a flow chart in Fig-<ref>. We start with a finite set of velocity realizations and associated probability weights (with a detailed discussion of sample choices in appendix A). Step 2: For each sample realization, the corresponding PF (or Markov) matrix, 𝐏_i, is constructed from the velocity field 𝐔_i for a specified time horizon (τ). Step 3: The set of contaminant tracking matrices, 𝐐_τ,i are constructed. Step 4: Operational scenarios are accounted for by applying sensor accuracy thresholding and constraints. This results in the set of matrices, 𝐐^*_τ,i Step 5: Each state is weighted by its volume to get the volumetrically scaled contaminant tracking matrices,𝐐^**_i . Step 6: The volumetric coverage, 𝐯̅_̅i̅, of each matrix is computed. Step 6: The expected coverage 𝔼[𝕍] is computed and the sensor locations calculated by finding the index of the maximum entry in the expected vector. 
 




§ RESULTS AND DISCUSSIONS

This section presents results of the algorithm for sensor placement under uncertain operating conditions for a set of benchmark buildings.


 §.§ Problem geometry and uncertain boundary conditions



We consider the IEA-Annex20 <cit.> benchmark problem in 2D and a furnished room in 3D. Fig <ref> shows the 2D geometry and the 3D complex building and the associated computational mesh. These meshes were chosen after performance of rigorous mesh convergence analysis. 2D and 3D simulations were performed on meshes with 9462 and 0.6M cells, respectively. 

We consider a scenario where the boundary condition (i.e. wall temperature) for the right wall (for both 2D and 3D) is uncertain. The right wall temperature is considered to be a random variable, ξ. We construct the distribution using a data-driven, location-specific approach as detailed in Sharma et. al <cit.>. We use the TMY weather data for Des Moines, IA, U.S.A and perform a whole year simulation for an identical domain in the building simulation software EnergyPlus <cit.>. We extract the right wall temperature from the output logs of the simulation. These are realizations of the random variable, ξ. We then utilize standard techniques – kernel density estimation <cit.> – to fit a probability distribution to the data. Fig. <ref> shows the data-driven probability distribution of the right wall temperature, as well as the associated cumulative distribution function and the inverse cumulative distribution function.







 §.§ Selecting the samples and checking for convergence in number of samples

The next step is to construct the flow realizations by selecting samples of the random variable, ξ (right wall temperature). A finite number of samples (and their associated probability) of the random variable are chosen so that they faithfully represent the true distribution of the random variable. The sampling of ξ is carried according to the procedure outlined in appendix-<ref>. We choose the number of samples to be M = 7. Since our final metric is the expected coverage, we perform a rigorous convergence analysis by increasing the number of samples, M, and evaluate the expected coverage for each M, until there is insignificant change in result (≤ 4%). The results of this convergence analysis is shown in Table <ref>, and indicates that M = 7 samples produces a sufficiently accurate representation of the random variable. While this convergence analysis was performed for the 2D case, our a posteriori analysis of the 3D results also indicated that this was a sufficient number of samples for the 3D case.





 §.§ Simulation of the flow fields and construction of the PF Operators
 
We perform non-isothermal simulations in 2D and 3D for each of the M = 7 cases identified. We utilize the OpenFOAM <cit.> simulation tool, specifically the buoyantboussinesqsimpleFoam solver. We set a convergence threshold of 10^-5 to solve the coupled Navier Strokes and energy equations. All walls have no-slip boundary conditions, with velocity inlet and zero-gradient outlet conditions for both 2D and 3D problem. All the walls except the right wall are considered well insulated. The inlet air temperature and the right wall temperature for each of the M = 7 sampling points are listed in Table-<ref>. The inlet air temperature is chosen such that indoor thermal comfort is maintained <cit.>. 


The computational mesh used for the numerical simulation significantly affects the results and hence grid convergence study is carried for both 2D and 3D geometries (please see appendix. <ref> for detailed mesh convergence results). The study resulted into the choice of 9462 cells in 2D and 0.6 Million cells for the 3D geometry. Following the computation of the velocity fields, we next turn to the construction of the PF matrix. As discussed in section-<ref> the size of the PF matrix is depend on the number of states/cells used in computation of flow field. We choose not to use the fine scale mesh used for CFD simulations to construct the PF matrix for each realization[Using the CFD mesh would result in a PF matrix of size 0.6M× 0.6M. Such a large matrix will require a lot of memory to store, and subsequent analysis will be infeasible]. We project the computed flow fields from the CFD mesh onto a coarser mesh containing around 70K cells. Therefore, the PF operator constructed is of size 70K×70K. A solver makeMarkovMatrixSteadyDiff is developed in OpenFOAM to construct the PF matrix by using the flow field obtained after solving NS equations. We make this solver freely available at <git@bitbucket.org:baskargroup/markovmatrixgeneration.git>.



Figure-<ref> shows the 2D velocity field contours computed by solving the steady NS equations for the sampled seven boundary conditions. It can  seen that various flow conditions are observed in the building by simply changing the right wall boundary conditions. The PF operator for each of the M = 7 cases is used in the sensor placement for the 2D problem. Figure-<ref> shows contours of flow field and temperature on the mid-plane cross-sections along with the sparse matrix representation of the PF operator. These PF operators are used to find the sensor locations for the 3D geometry on the expectations based sensor placement framework developed in section-<ref>.




 §.§ Sensor placement accounting for uncertainty: 2D case



  §.§.§ Validation of the contaminant transport via PF operator
2D

To validate the constructed PF operators we compare the contaminant transport, Φ, computed by the PDE Eq. (<ref>) and Markov matrix Eq. (<ref>). While we validate each of the M = 7 PF matrices against a PDE based result, we show results for a single case for brevity. Figure-<ref>(a) shows the initialized scalar map where red represents a scalar concentration of one and blue corresponds to zero. The scalar transport PDE equation is solved to a time horizon of 50 sec using the OpenFoam solver, and compared with the Markov approach. Figure <ref>(a-b) compares the scalar concentration contours after 50 sec by PDE and Markov approach. It can be seen that they closely match each other. The concentration profiles are also plotted  and compared  about the mid-planes of the building which is shown in fig.<ref>. The Markov results are represented as dots, and very closely match the PDE based predictions with an L_2 error of 10^-4. The validation shows the accuracy of the Markov approach and effectiveness of the matrix-vector product based approach in predicting the transport of scalar concentration.







  §.§.§ Sensor locations with no-constraints

We identify sensor locations under uncertain conditions for the no constraint scenario. We place four sensors in the domain with a sensing time of t= 60 sec. The accuracy of the sensor used for thresholding the contaminant matrix is ϵ_acc = 0.01%. The optimal sensor location are shown in Fig. <ref>(a). The probability coverage for each sensor is shown in Fig.<ref>(b-e). 
For the sensors the associated volumetric coverage for the sampled seven realization and the extent of cover is shown in Fig. <ref>(a-b). The total expected coverage for all the realization using four sensors is  83.65%. Having illustrated the approach for a 2D benchmark, we next deploy this approach for a complex 3D scenario. 






 §.§ Sensor placement accounting for uncertainty for a complex building: 3D case

Next we discuss the results of sensor placements for the 3D complex furnished building geometry. The constructed PF operators and the flow fields are already shown in section-<ref>. We first validate the accuracy of the PF operator by comparing PF based contaminant transport predictions with PDE based predictions.



  §.§.§ Validation of PF based contaminant transport

We perform validation similar to the 2D PF operator validation case discussed in section-<ref>. The transport of contaminant, Φ, is performed using both PDE solutions and PF operator. The results are shown for one of the constructed PF operator but all the constructed PF operators are validated in similar manner before using them for sensor placement. Figure-<ref>(a) shows the initial distribution of the contaminant at t=0 where the contaminant density of one is marked with red and zero where the contaminant is absent. Figure-<ref>(b) shows iso-surface of the evolution of the contaminant from the solution of PDE in the domain at t=50sec. Figure-<ref>(c) show the iso-contours of the contaminant by the PF operator at t=50 sec. It can be observed that two iso-contours are indistinguishable. To further validate the PF operator, we monitor the concentration profile along the three center axes passing, as shown in Figure-<ref>. It can be seen that they exactly overlap each other which shows the accuracy of the PF operator in transporting the contaminant. 







  §.§.§ Sensor placement with no constraints

We first place a fixed number (k = 4) of sensors to maximize the expected coverage without any sensing and location constraints. Fig.<ref>(a) shows the sensor location under this case. It can be seen that to maximize the expected coverage the first and third sensors are placed inside the occupied region. Another interesting observation is related to the placement of the first sensor. For the deterministic case, the first sensor is usually placed close to the outlet, as has been extensively reported by <cit.>. Placing the first sensor close to the outlet maximizes the observability of the domain, as all the air (and contaminant) must exit from the outlet. In contrast to the deterministic case, maximization of the expected coverage is not achieved by placing the first sensor near the outlet but (for this geometry) is inside the occupied zone of the building. 





  §.§.§ Sensor placement under location constraints

The location constraints are particularly attractive for critical applications such as detection of CBW (chemical and biological weapons) and TID (transmission of infectious diseases). As most of the CBW and TID sensors are expensive and can be damaged if placed in the occupied space, this constraint forces the sensors to be placed outside the occupied zone. Fig.<ref>(b) shows the obtained sensor location for four sensors. 



  §.§.§ Sensor placement with sensing and location constraints

In case of indoor air quality applications, the sensing of the contaminant concentrations
in the occupied zone is much important than outside the occupied zone. Therefore a sensing constraint is added while determining the optimal location of the sensor. Similar to the case of CBW and TID scenarios, it is important that the sensor must not be placed where the occupant can affect or damage it. For this reason a location constraint is also incorporated so that the sensor must be placed outside the occupied zone. Fig.<ref>(c) shows the sensor location in the isometric view and top view. In this case sensors tend towards to be placed near the occupied zone or at the outlet. 



  §.§.§ Discussion of quantitative measures

 We finally compute quantitative measures for the placed sensors. Under the no-constraint scenario, the optimal placement of sensors results in an expected coverage of 31.4 % of the total volume. That is, under all possible flow condition scenarios produced by an uncertain wall boundary condition, the sensors are expected to cover 31.43 % of the domain. Interestingly, under location constraints, the expected coverage does not change much from the no-constraint case. This is primarily because of the increase in coverage sensor-2 provides for realization-2. With the addition of both location and sensing constraints there is a  decrease in the expected coverage to 15.0 % of the total volume. However, since in this case we are only interested in the occupied zone, the sensor placed with sensing and location constraints cover 21.50 % of the occupied space. 

For each of the M = 7 realizations, we compute the volumetric coverage for each realization associated with each sensor. Figure-<ref> shows plots of these quantitative measures, where fig-<ref>(a-b) corresponds to the sensor placement with no-constraint on placement. (c-d) corresponds to the constraints on the location as well as sensing and (e-f) corresponds to the sensor placement with a only constraint on the location of the sensor. By comparing the extent of cover for fig-<ref>(a-c-d) it can be seen that the sensor positions obtained under no constraint result in high extent of cover for all realizations. While for the cases with constraint on placement the extent of cover varies significantly with each realization. Comparing the volume cover for fig-<ref>(e-d-f) it can be observed that under no constraint sensor placement the first sensor performs majority of sensing for each realization, while in case of sensing and location constraint the sensing is well distributed among the four sensors.   





§ CONCLUSION

The work introduces a method for accounting for weather uncertainty during optimal sensor placement. The method is based on constructing the Perron-Frobbenious operator from the flow field using a set based approach. We use the theoretical background of uncertainty quantification to formulate the sensor placement problem. The weather affecting the flow field is accounted as uncertainty in the current work. We utilize a data-driven approach to construct the uncertain boundary condition as a probability distribution. We then select a finite number of representative samples from this distribution. The samples and their associated probability of occurrence are chosen such that the expected value of these samples is identical to the full distribution. Full CFD simulations are performed corresponding to this finite set of boundary conditions. The velocity fields are then used to construct individual PF operator which is used in the algorithm to identify the sensor layout for the domain.

This approach can naturally account for any uncertainty to determine the optimal sensor layout. By using a dynamical systems approach we gain two key advantages: (a) we circumvent the numerical complexity and computational challenges with complementary optimization/inverse based approached. Additionally, incorporating effects of uncertainties is non-trivial in other methods, (b) the PF based approach produces a mathematically rigorous and globally optimal sensor layout. This is in contrast to other optimization based approaches which can only identify local optima. We illustrate the approach for problems in both 2D and 3D building geometry. We envision that the developed framework can be used for critical applications including sensor layout for CBW scenarios, transmission of infectious diseases, and indoor air quality. The contributions from this paper also provide ability for designers, engineers, and researchers to better understand the probabilistic coverage maps for each sensor location. 



§ ACKNOWLEDGEMENT

UV acknowledges partial support from NSF CAREER 1150405. BG acknowledges partial support from NSF CAREER 1149365.


  


§ APPENDIX



 §.§ Sampling points and associated weights

We use sampling theory to choose a finite set of realizations of a random variable, ξ, with a probability distribution, ρ(ξ).  Given an output variable, 𝐯(ξ), that depends on the random variable, the expected value of this output variable is given by

    𝔼[𝐯] = ∫𝐯(ξ) ρ(ξ) dξ

Now, the finite number of realizations, {ξ_1, ⋯, ξ_M } are chosen such that they can accurately recreate the expected value. That is

    ∫𝐯(ξ) ρ(ξ) dξ = ∑_i=1^M𝐯(ξ_i)θ_i

where 𝐯(ξ)_𝐢 are the outputs evaluated at the sampling points ξ_i, and θ_i are the probabilistic weights (defined below). We use the approach detailed in Zabaras and Ganapathysubramanian <cit.> to sample from a distribution[In a data-driven context, we do not necessarily have a parametrized distribution for ξ. The usual approach is to use a kernel density estimator to construct a non-parametric representation of the distribution from data.]. Essentially, we sample points uniformly from the inverse cumulative distribution function (iCDF) of ξ and increase the number of points until the expected value has converged. This ensures that we have chosen representative points which provides a statistical representation of the actual distribution of the random variable. In case of more number of random variable, sampling techniques such as latin-hypercube sampling (LHS),importance sampling and sparse grid collocation can be used <cit.>.       

The determination of the weights θ_i is performed using interpolation theory. The output variable, 𝐯, which is a function of ξ is represented as a sum of basis functions N_i(ξ) (defined over the support of ξ). That is

    𝐯(ξ) ≈∑^M_i=1𝐍_i(ξ)𝐯_i

where 𝐯_i are the outputs evaluated at the points ξ_i. 
Inserting this representation into Eq. <ref> results in

    𝔼[𝕍]    = ∫∑^M_i=1𝐍_i(ξ)𝐯_iρ(ξ) dξ
       =  ∑^M_i=1𝐯_i∫𝐍_i(ξ) ρ(ξ) dξ
       = ∑^M_i=1𝐯_iθ_i

Thus, θ_i = ∫𝐍_i(ξ) ρ(ξ) dξ. 

To illustrate the approach we consider ξ to be a standard normal variable exhibiting a Gaussian distribution (𝒩(0,1),μ = 0.5,σ=0.05). We evaluate θ_i and compare the true integral (Eq. <ref>) with the approximate evaluations using M = 7 sampling points, for three different functional forms of 𝐯. Linear basis functions were used to approximate 𝐯(ξ) in Eq.(<ref>). Table-<ref> show the comparisons with three different functions. It can be seen that the approximation integral closely matches the exact integral up to two decimal places. 




§ GRID CONVERGENCE

 The grid convergence study is performed  using cdf-0.1 boundary conditions shown in table-<ref> for the 3D geometry. The temperature profiles are shown in Fig.<ref> plotted along the x-axis passing through the room center and along an arbitrary chosen y-axis close to cabinet ((7.5,0,0.5),(7.5,3.0,0.5)) for the 3D geometry(fig. <ref>). It can be seen that for all the meshes the X-centerline profile closely overlap each other. While in case of Y-line the difference are merely ± 0.2 K between the fine 4.0M mesh and the chosen 0.6M mesh. Therefore, all the CFD  computations to compute the flow field were carried using 0.6M mesh.




model1-num-names




Usteady/unsteady flow field
𝐏Markov Matrix
Ωcomputational domain
Φcontinuous form scalar contaminant
S_Φsource of contaminant

Ŝ_t_i,t_i+1 discrete form of source location
ttime [sec]
t_ffinal time of simulation [sec]
δ ttime associated with Markov matrix
τsensor sampling time [sec]
𝐐_τcontaminant tracking matrix
ϵ_accsensor accuracy provided by manufacturer
𝐐^*_τthreshold contaminant tracking matrix

𝔼expectation operator
ξrandom variable


ρprobability density function
ρ̂approximated probability density function
Mnumber of realizations for random variable ξ
Θthe weights associated with each sample

