**4. Simulation Design**

A model of the UBO was first created with the aid of SketchUp [34], a 3D modeling software. Using dimensions taken from the building, the single story layout was replicated. Utilizing the plug-in from OpenStudio, thermal zones, boundaries, and interactions were defined. This model was then exported as an input file (IDF) for EnergyPlus [35], an open-source energy simulation program that has been developed by the Building Technologies Office (BTO) under the U.S. Department of Energy (DOE). Using EnergyPlus, HVAC equipment was added to the model based on the equipment present in the UBO. While EnergyPlus excels at modeling and simulation, it is not immediately accessible for controller development and implementation. As such, controllers for the AHU and the zone level VAVs were created in MATLAB. To enable co-simulation between EnergyPlus and MATLAB, two programs were used. The first was MLE+ [36], an open-source MATLAB toolbox for creating the necessary configuration files and providing functions to connect EnergyPlus and MATLAB with an easy-to-use graphical interface. MLE+ utilizes the Building Controls Virtual Test Bed (BCVTB) [37] as the communication backend between EnergyPlus and MATLAB, providing the co-simulation functionality.

The individual zones of the model were setup as defined in the building layout and exterior doors and windows were placed as accurately as possible. The AHU and zone VAVs were added using EnergyPlus' HVACTemplate objects. A central electric chiller was also added to supply chilled water to the AHU. The chilled water output temperature is regulated to 45 °F, the same as the supply chilled water temperature for the UBO.

Currently, EnergyPlus does not offer modeling of pressure with variable air volume systems, thus another solution was necessary to simulate the UBO's control and physical limitations of end static pressure and flow in the AHU. To accomplish this, the authors assumed that the dynamics of the fan speed and end static pressure were fast enough compared to the simulation timestep (1 min) to be considered instantaneous. Additionally, the assumption that the fan would supply the requested end static pressure, constrained by the physical limitations of the AHU and ducting, was made. To determine this constraint, data from the real UBO building were analyzed and a maximum possible work performed by the fan was calculated (910 W). During the optimizations, the constraint is calculated by using Equation (19).

To determine the total volume flow through the AHU, individual models of the zone VAVs were generated from data based on VAV damper position and end static pressure. The effect of outdoor air temperature was also considered on the VAVs as the AHU draws in outdoor air, but was shown to be minimal. This is most likely due to the fact that AHU is able to meet its discharge air temperature setpoint, even at varying flows, effectively isolating the VAV supply air from the outdoor air conditions. The flows from the VAV models are summed to obtain the total flow through the AHU, assuming minimal duct losses. The constraint can be written as:

$$910 \ge 0.1175 \cdot \frac{q\_{AHII} \cdot P\_{EDS}}{\mu\_f \cdot \mu\_b \cdot \mu\_m} \tag{19}$$

Thus, the optimization will only ever choose an end static pressure setpoint that is physically achievable by the system. This end static pressure is then passed to the VAV models which, along with the commanded damper positions, produce individual zone air volume flows.

The overall control hierarchy can be seen in Figure 9. The supervisory controller supplies the zone temperature setpoints ( *<sup>T</sup>*<sup>∗</sup>*ZONES*) to the respective PID reference inputs and the end static pressure setpoint ( *<sup>P</sup>*<sup>∗</sup>*EDS*) to the VAV models. The discharge air temperature setpoint ( *<sup>T</sup>*<sup>∗</sup>*AHU*) is supplied directly to EnergyPlus as the control for the chilled water valve is implemented using an appropriate setpoint manager within EnergyPlus. The zone temperature error is fed to the cooling PID controller which outputs a desired percentage of maximum flow setting (0–100%). If the minimum flow setting for a zone is zero, then the signal remains unchanged. The desired flow percentage is then passed to the flow PID controller which then produces a desired damper position. This damper position is used in the VAV models as previously described. The outputs from the VAV models, desired zone air volume flows, are converted to air flow fractions and then sent to zone VAVs in EnergyPlus where the room dynamics are simulated for one timestep.

**Figure 9.** Control hierarchy used in simulation of the UBO.

For the supervisory controllers, a prediction model is necessary to determine the future zone temperatures as the setpoints are optimized. A previously developed modeling algorithm was employed to generate models for the individual zones [38]. One of the main reasons this algorithm was chosen is because of its ease in producing models, requiring the user to only select input and output data. The full details of this process are beyond the scope of this publication, but more information can be found in [38]. Briefly, data from the EnergyPlus simulation are analyzed to develop discrete, linear models of the selected parameters. These models take the form of:

$$\begin{aligned} \mathbf{x}\_i(k+1) &= A\_i \mathbf{x}\_i(k) + B\_i \mathbf{u}\_i(k) + K\_i e\_i(k) \\ \mathbf{y}\_i(k) &= \mathbf{C}\_i \mathbf{x}\_i(k) \end{aligned} \tag{20}$$

where *A*, *B*, *C*, and *K* are the identified system matrices, *u* is the model input vector, *y* is the predicted output, and *e* is the error defined as the difference between the current value and the previous predicted value. The algorithm automatically identifies significant coupling interactions between zones and includes the respective zone temperatures (*Tzone dist*.,*j*) as measured disturbances. In addition to the temperature of these disturbance zones, other parameters such as outdoor air temperature (*TOA*), outdoor air relative humidity (*RhOA*), AHU discharge air temperature (*TAHU*), end static pressure (*PEDS*), and zone temperature setpoints (*T*<sup>∗</sup>*zone*,*<sup>i</sup>*) are used as inputs to ARX, ARMAX, Output Error, and Box-Jenkins modeling methods with the output being the respective zone temperatures (*Tzone*,*<sup>i</sup>*). The best fitting model is selected and then the individual models are combined into a centralized model of the entire system. Steady-state predictions from this centralized model are then used by the supervisory controller to optimize the UBO's 12 setpoints (discharge air temperature setpoint, end static pressure setpoint, and 10 zone temperature setpoints) to minimize the economic cost functions previously described. The optimization occurs every 15 min with the *R*, *S*, and *T* matrices updating based on current operating conditions to reflect the time varying nature of the system dynamics.

Steady-state relationships were chosen over dynamic relationships for the initial implementation and validation of the proposed economic objective function strategy. This choice takes advantage of the fact that in the building energy systems, the control variable's dynamics (AHU discharge air temperature, end static pressure, VAV damper position) change quickly versus the other system variables (outdoor air temperature, outdoor relative humidity, room temperature), which change

relatively slowly over the optimized timestep. The exact model inputs and outputs are shown in Table 3.


**Table 3.** Inputs and outputs for the generated zone models.

Models were identified for two operational cases: (1) where the zone VAV damper still had actuator range (i.e. the damper was not fully open); and (2) where the zone VAV damper is fully open. The separate models were necessary as the effect of the model inputs varies greatly between the two cases. In Case 1, the effect of *<sup>T</sup>*<sup>∗</sup>*AHU* and *<sup>P</sup>*<sup>∗</sup>*EDS* are minimal compared to *<sup>T</sup>*<sup>∗</sup>*zone*,*i*. This is due to the fact that the VAV damper is actuated by a PID controller with *<sup>T</sup>*<sup>∗</sup>*zone*,*<sup>i</sup>* as the reference. The PID controller is able to reject changes in *<sup>T</sup>*<sup>∗</sup>*AHU* and *<sup>P</sup>*<sup>∗</sup>*EDS* by changing the damper position to achieve *<sup>T</sup>*<sup>∗</sup>*zone*,*i*. However, in Case 2, when the damper is fully open, *<sup>T</sup>*<sup>∗</sup>*AHU* and *<sup>P</sup>*<sup>∗</sup>*EDS* become the inputs of significance as they directly effect the zone temperature, determining the amount and the temperature of the incoming conditioned air.

A decision process was necessary to determine when each model should be used. This decision was made based off of two conditions. The first determined if the predicted zone temperature using the Case 1 models was greater than the prediction from the Case 2 models. This condition served to verify whether the predicted temperature from the Case 1 models was currently achievable with the state of the AHU. If the Case 1 predicted temperature was lower than the Case 2 predicted temperature, then the VAV would not be able to achieve the Case 1 predicted temperature with the current *<sup>T</sup>*<sup>∗</sup>*AHU* and *<sup>P</sup>*<sup>∗</sup>*EDS* values. The second checked if the current VAV damper position was less than 95%, or, in other words, if the VAV still had actuator range of the damper. The value of 95% was used as opposed to 100% to serve as a threshold and help prevent the system from oscillating between cases. If both these conditions were true, then the Case 1 models were used; otherwise, the Case 2 models were used. To help ensure smooth transfer between the two models and more accurate predictions, errors were calculated between the previous predictions the measured temperatures and included in the current prediction.
