Development of a variable refrigerant flow system emulator to host the second World Championship in Cybernetic Building Optimization

Abstract

This paper introduces an emulator developed to assess the operational characteristics of variable refrigerant flow (VRF) for the second World Championship in Cybernetic Building Optimization (WCCBO). The emulator utilizes BACnet communication and provides a realistic simulation for evaluating energy conservation and comfort in buildings with VRF systems. Herein, we discuss the building specifications for calculation, operation of the emulator, calculation method for the vertical temperature distribution, dissatisfactory rate calculation, and consideration of championship rules. Notably, a thermal preference model was developed for predicting the thermal sensation of occupants and adjusting the dissatisfaction rate based on their ability to control the indoor environment. In addition, we developed a method for calculating the total scores during the WCCBO and devised a suitable network structure for score registration. The results of this study highlight the importance of personal control and impact of drafts on thermal comfort. This paper presents a meta-level methodology to determine the best operation for each specific VRF system for a future world championship, whose planning is underway and will be conducted to accumulate information on optimal operations.

KEYWORDS:

• Emulator
• VRF
• HVAC
• thermal comfort
• occupant

Introduction

This paper reports the development of an emulator to assess the skill of operating variable refrigerant flow (VRF) systems. The authors developed a similar emulator for central heat source systems (Togashi and Miyata 2019) and reported the results of a championship on optimization using an emulator (Togashi, Miyata, and Yamamoto 2020). A second championship is currently being planned for the VRF system, and the emulator reported in this study was developed for use in this championship. The amount of energy consumed by buildings is significant (DOE 2015; European Parliament 2010; METI 2015), and many researchers have reported the importance of optimizing the operation of heating, ventilation, and air conditioning (HVAC) systems to reduce energy consumption (Chau, Leung, and Ng 2015; de Wilde 2014; Khoury, Alameddine, and Hollmuller 2017; Sartori and Hestnes 2007; Verhelst et al. 2017). To improve energy efficiency of buildings, a device that can correctly evaluate the skill of HVAC system operation is required. Thus, we propose the application of a highly realistic simulator (hereinafter referred to as ‘emulator’) that returns a response that is almost identical to that of a real building. The history of the emulator concept, its application for evaluating the skill of HVAC system operation, and the functions required for emulators have already been fully discussed by Togashi and Miyata (2019) and will not be repeated here.

Notably, HVAC may be controlled manually via humans or through computer algorithms; however, in this study, regardless of who controls it, the ‘skill of HVAC system operation’ refers to whether the operation is superior in two respects: energy performance and comfort.

The novelty of this study is that it focuses on VRF systems. As reported in detail in Togashi and Miyata (2019), the development of emulators has been extensively conducted (Blum et al. 2021; Bushby et al. 2001; 2010; IEA 1996; 1997; Lebrun and Wang 1993; Vaezi-Nejad et al. 1991); however, till now, thermal environment systems, including VRF systems, have not been emulated. Nevertheless, in recent years, the percentage of buildings with VRF systems has continually increased (Aynur 2010; Lin et al. 2015; Zhang et al. 2019). Education of engineers is one effective use of emulators, and specific examples of education through simulation have been reported by Neuman (2011; 2012), Serra et al. (2017), Gentile, Kanters, and Davidsson (2020), and Madina, Pratiwi, and Tundono (2021). It would be of great significance to develop an emulator system that is useful for education on the operation of VRF systems, which are being widely adopted in society.

The authors held a championship in 2019 to compete in building optimization using the developed emulator (Togashi, Miyata, and Yamamoto 2020). In this paper, this championship is referred to as the first WCCBO (The first World Championship in Cybernetic Building Optimization (WCCBO). In the first WCCBO, 34 participants repeated 339 annual virtual operations using emulators to compete in terms of energy conservation and comfort. The buildings evaluated in the first WCCBO had central heat source systems. Therefore, we plan to organize a new championship for buildings with VRF systems. Hereafter, this will be called the second WCCBO (The Second World Championship in Cybernetic Building Optimization (WCCBO). The VRF system emulator reported in this paper was specifically designed for use in this championship.

Figure 1 shows the overall schematic of this study. The upper panel shows the stand-alone VRF System Emulator, and the lower panel shows the 2nd WCCBO to be held using the emulator.

Figure 1. Schematic of the study objective.

In this study, we defined the ‘VRF System Emulator’ as a composite model that includes the VRF system model and the models of the occupant and the building (Figure 1, upper panel). Individually, these three components have been extensively studied. However, for use as an emulator, two features are required. One is the requirement of a generic interface to receive external control. The other is to the ability to model the interaction effects of components and express tradeoffs between energy performance and comfort. The implementation of these features is a major feature of this study and is described below, along with the role of each Section of this paper.

To increase the versatility of the emulator, BACnet communication is used to operate the VRF system model by the HVAC Operator. The significance of such an emulator operation using BACnet communication has already been reported in a previous study (Togashi and Miyata 2019). The development and accuracy verification of the physical model of the VRF system has also been reported (Togashi and Satoh 2021). However, the type of BACnet signals (or BACnet objects) that the VRF system should accept is yet to be investigated. This is discussed in Section 3.

When the HVAC Operator gives a command via BACnet to the VRF system model, the VRF system and Occupant models output the energy consumption and the dissatisfied rate, respectively. These outputs change according to the HVAC Operator’s commands. The smaller the value of both, the better the skill of HVAC system operation.

The interactive effects of the VRF system, occupant, and building must be represented as per the real situation to assess whether the HVAC operation is well or poorly done. This is because, typically, there is a tradeoff between energy performance and comfort. For example, reducing ventilation improves energy performance because reduces the outside air load; however, it reduces comfort as it reduces air cleanliness. In another example, reducing airflow such that the office worker does not feel a draft may reduce the heat exchange efficiency inside the indoor unit. To correctly represent such trade-offs, calculations are not included in the general building (heat load calculation) or occupant models must be introduced. They are as follows.

The building model must be able to calculate the change in the airflow near the office workers and the temperature distribution created owing to the air blown by the indoor unit. However, such airflow is generally not output in heat load calculation simulation models. Considering the purpose of the emulator, it is of course assumed that time-consuming methods such as computational fluid dynamics (CFD) cannot be used. A solution to this problem is proposed in Section 4.

The occupant model must be able to model the thermal dissatisfaction and the dissatisfaction about elements that have tradeoffs with HVAC energy consumption, i.e. airflow, temperature distribution, cleanliness. It also needs a method to integrate these dissatisfactions. These methods is reported in Section 5.

In many real-world VRF systems, the occupant uses a remote controller to change the settings of the indoor unit. Such individual operability affects the psychology of the occupant as well as energy and comfort. Methods for modeling such behaviour of the occupant and evaluating its psychological impact are reported in Section 5.

The final goal of this study was to propose a method to determine the HVAC Operator who exhibited the best skill of HVAC system operation using the VRF System Emulator developed by employing the above method. This method is shown in Figure 1 lower panel, where the VRF System Emulator is distributed to each of the n HVAC Operators. In this figure, the HVAC Operator is referred to as the Competitor.

Each Competitor operated the VRF system model in a different way, thus outputting different energy consumption and dissatisfied rate. The network structure for collecting the results of these calculations while preventing tampering is described in Section 6. The ranking of the results of the two calculations: energy consumption and dissatisfied rate, is a major point of discussion. If there is only one indicator, they can be simply ranked in order of decreasing value. As the two indicator (energy performance and comfort) with different units exhibit a tradeoff relationship, they must be integrated in a fair manner. The calculation method for this integrated indicator is discussed in Section 6.

As mentioned above, the VRF System Emulator is a complex program with many factors influencing each other; thus, validating the response to any HVAC operation is challenging. Therefore in Section 7, we test several typical HVAC operations and interpret the results.

Specifications of building and HVAC system to be emulated

For the calculations, we used the same building as that used in the first WCCBO. The corresponding floor plan and air-conditioned zone are shown in Figure 2. The office building was used as a reference for the Japanese Energy Conservation Law (see SHASE (2016) and Ono, Ito, and Yoshida (2017) for details on the specifications). In the first WCCBO, a central heating source system was operated, and thus, the entire seven-story building was included in the calculations. By contrast, in the second WCCBO, a VRF system will be operated, and therefore, only half of one floor was included in the calculations in this study. Further, contrary to the first WCCBO, in which an evaluation period of one year was set because of the extremely large thermal storage tanks, in the second WCCBO, the evaluation period will span one week each in summer and winter because such a large thermal lag is not expected in the case of VRF systems.

Figure 2. Floor plan and HVAC zoning of the simulation target.

Because the effects of solar radiation vary with the azimuth, the optimal method of HVAC operation also varies with the azimuth. Thus, VRF systems in two rooms, each on the north and south sides, were targeted for operation. These rooms are occupied by a different tenant. The behaviour of each office worker was calculated independently by entering and leaving the building at different times. The details of this office worker behaviour model are the same as in the first WCCBO, please refer to Togashi and Miyata (2019) for details on the calculation method.

As shown in Figure 2, the room is divided into nine air-conditioned zones, each with its own indoor units. The five perimeter zones and four interior zones have separate outdoor units; therefore, there are two VRF systems in the north tenant and two in the south tenant, for a total of four VRF systems. The four colours in the figure represent four different VRF systems.

In addition to VRF indoor units, heat recovery ventilation units are installed in each zone. These units could bypass the outdoor air without heat exchange. The air blown out by the ventilation units are directly supplied to the zones.

Humidifiers with a saturation efficiency of 95% and an airflow rate of 9.4 CMH per floor area are installed in each zone, and they start and stop to maintain the relative humidity in the zone between 40% and 50% during winters. Because they are drip-type humidifiers, heat equivalent to the latent heat from humidification is added to the sensible heat load of the indoor units.

Table 1 and 2 list the specifications of the outdoor and indoor units of the VRF system, respectively. Table 3 lists the specifications of the ventilation units. Table 4 lists the indoor units and ventilation units in each zone. JIS in the table refers to the rated conditions for disclosing performance as defined by the Japan Industrial Standard. Togashi and Satoh (2021) presents the details related to the initialization of the VRF system model from these specifications.

Table 1. Outdoor unit specifications.

-	Conditions		V1	V2	V3	V4
Cooling	Capacity	JIS nominal	40.0	22.4	33.5	22.4
	[kW]	JIS mid-load	18.0	10.1	15.1	10.1
		JIS mid-load and mid-temperature	18.8	10.6	15.7	10.6
	Electricity	JIS nominal	12.5	6.07	9.74	6.07
	[kW]	JIS mid-load	3.50	1.89	2.59	1.89
		JIS mid-load and mid-temperature	2.74	1.55	2.12	1.55
Heating	Capacity	JIS nominal	45.0	25.0	37.5	25.0
	[kW]	JIS mid-load	20.3	11.3	16.9	11.3
	Electricity	JIS nominal	13.1	6.32	10.0	6.32
	[kW]	JIS mid-load	3.89	2.06	2.94	2.06
Fan	Air flow rate [m^3/min]		210	218	187	218
	Electricity [kW]		0.58	0.52	0.42	0.52
Pipe	Nominal length [m]		7.5	7.5	7.5	7.5
	Comparison length [m]		100	100	100	100
	Correction factor [-]		0.885	0.880	0.895	0.880

Table 2. Indoor unit specifications.

Indoor unit type	C56	C71
Nominal cooling capacity [kW]	5.6	7.1
Nominal heating capacity [kW]	6.3	8.0
Air flow rate [m^3/min]	15.5	22.0
Electricity [kW]	0.043	0.072

Table 3. Heat recovery ventilation unit specifications.

Ventilation unit type	HX150	HX250
Ventilation air flow [CMH]	150	250
Sensible efficiency [%]	74.0	70.0
Latent efficiency [%]	64.5	63.0
Electricity [kW]	0.075	0.100

Table 4. Type of indoor-units and ventilation units in each zone.

Zone name	N1	N2	N3	N4	N5	N6	N7	N8	N9
I/U name	V3-1	V3-2	V3-3	V3-4	V3-5	V4-1	V4-2	V4-3	V4-4
I/U type	C71	C56	C56	C71	C71	C56	C56	C56	C56
HEX type	HX150	HX150	HX150	HX250	HX250	HX250	HX250	HX250	HX250
Zone name	S1	S2	S3	S4	S5	S6	S7	S8	S9
I/U name	V1-1	V1-2	V1-3	V1-4	V1-5	V2-1	V2-2	V2-3	V2-4
I/U type	C71	C71	C71	C71	C71	C56	C56	C56	C56
HEX type	HX150	HX150	HX150	HX250	HX250	HX250	HX250	HX250	HX250

The capacity of the indoor units was selected based on the results of the maximum heat load calculations. The cooling capacity was approximately 200 W/m² per floor area. This capacity is appropriate because the average cooling capacity of VRF systems in office buildings of 300–2000 m² in Tokyo in FY2022 is 235.4 W/m² according to a survey by the Japanese Ministry of Land, Infrastructure, Transport and Tourism (Miyata and Miki 2024). However, in recent years, the amount of heat generated in offices has been decreasing owing to the change to LED lighting and use of more efficient computers. Therefore, a cooling capacity of 200 W/m² is excessive for the actual heat load. In reality, many VRF systems are operated at a lower partial load. The operation of the VRF system under such conditions is another issue that we wish to answer using this emulator. Note that the internal heat generation assumed in the emulator’s building model was 7.1 W/m² for lighting, 40 W/person per occupant staying in office for plugs, and 3 W/m² for standby power.

VRF system model

Physical model

Many researchers have proposed simulation models for VRF systems (DOE 2019; Fujii et al. 2009; Hong et al. 2016; Matsumoto et al. 2015; Saito and Jeong 2012; Sato et al. 2008; Shinagawa et al. 2018; Torregrosa-Jaime et al. 2018; Zhou et al. 2008; Wang et al. 2023; Xiao et al. 2023). To use the model as an emulator and test skill of HVAC system operation, it should be capable of performing a simulation for time periods longer than a week at high speed and express the energy characteristics when the office worker operates each indoor unit individually. However, no model exhibits all these characteristics.

Thus, the authors proposed a model of a VRF system that can calculate annual energy in a short time and can reflect differences in the input conditions for each indoor unit (Togashi and Satoh 2021). Moreover, this study further developed this model to facilitate communication using BACnet protocol.

The BACnet protocol was used to fulfill two aims. The first was to allow the user to experiment with various controls. The user’s only constraint is to send commands via the BACnet, and the control content can be calculated using any method. The system can be operated manually using general-purpose BACNet communication software, or users can create their own complex control software. Because any programming language can be used, it is easy to apply machine learning techniques for control using Python. The second objective is to enhance the application scope of the control technology developed in the competition to real buildings. If the control system is specialized for a particular simulation software program, then it cannot be used directly to solve real-world problems, even if knowledge is accumulated through competition. If the emulator interface is similar to the control system of a real VRF system, then it can be directly used for the optimization of real buildings.

Accordingly, two issues are addressed in the following section: (i) identification of control items for the emulator, and (ii) formulation of a suitable method to synchronize time with that of the emulator, whose operation time is different from real time.

Determination of VRF monitoring parameters and control items

As mentioned previously, an emulator should have an interface similar to that of a real-world VRF system. Thus, the BACnet communication specifications published by the current major VRF manufacturers (Daikin 2011; Johnson Controls-Hitachi Air Conditioning 2017; Mitsubishi Electric 2019; Mitsubishi Heavy Industries 2019; Panasonic 2023; Toshiba Carrier 2017) were surveyed to identify the items that each company could monitor and control.

Table 5 presents the items monitored and controlled by each company. The rightmost column lists the items monitored and controlled by the emulator based on the survey results. Shaded items are provided by most manufacturers, and are therefore important. However, since the emulator is intended to evaluate energy and comfort, ‘Alarm’ were deemed unnecessary and removed. Although a dirty filter can affect energy consumption by reducing the heat exchange efficiency, the ‘filter sign signal’ was deleted because this championship was not intended for long-term operation, where the filter can become dirty.

Table 5. Monitor or control item supported by VRF manufacturers.

No.	Target			Manufacturer						EM
				D	ME	T	P	JH	MH
1	Indoor unit	On/Off	Control	x	x	x	x	x	x	x
2			Monitor	x	x	x	x	x	x	x
3		Operation mode	Control	x	x	x	x	x	x	x
4			Monitor	x	x	x	x	x	x	x
5		Room temperature setpoint	Control	x	x	x	x	x	x	x
6			Monitor	x	x	x	x	x	x	x
7		Room temperature	Monitor	x	x	x	x	x	x	x
8		Room humidity	Monitor							x
9		Alarm	Monitor	x	x	x	x	x	x
10		Malfunction code	Monitor	x			x	x	x
11		Connection error	Monitor				x	x	x
12		Filter sign signal reset	Control	x	x	x	x	x	x
13		Filter sign signal	Monitor	x	x	x	x	x	x
14		Fan speed or volume	Control	x	x	x	x	x	x	x
15			Monitor	x	x	x	x	x	x	x
16		Airflow direction	Control	x	x	x	x		x	x
17			Monitor	x	x	x	x		x	x
18		Indoor fan status	Monitor	x			x
19		Aux. heater status	Monitor	x
20		Ventilation On/Off	Control			x
21			Monitor			x
22		Ventilation mode	Control	x	x	x
23			Monitor	x	x	x
24		Ventilation amount	Control	x	x	x
25			Monitor	x	x	x
26		Forced Thermo-off	Control	x			x
27			Monitor	x			x
28		Thermo on/off status	Monitor				x
29		Total operating hours	Monitor			x
30		Indoor unit heat demand	Monitor			x
31		Trend log (operation mode)	Monitor			x
32		Trend log (temperature)	Monitor			x
33		Trend log (operating hours)	Monitor			x
34		Trend log (heat demand)	Monitor			x
35		Trend log (setpoint temperature)	Monitor			x
36		Fan electricity	Monitor							x
37		Heat load	Monitor							x
38	Outdoor unit	Compressor status	Monitor	x
39		Silent mode	Control					x
40			Monitor			x
41		Capacity control	Control					x
42			Monitor
43		Fan and compressor electricity (or gas)	Monitor	x					x	x
44		Heat load	Monitor							x
45	Whole system	Remote controller Permit / Prohibit	Control	x	x	x	x	x	x	x
46			Monitor	x	x	x	x	x	x	x
47		Energy saving mode	Control	x			x
48			Monitor	x			x
49		Forced system stop	Control	x					x
50		Forced evaporating/condensing temperature control	Control							x
51			Monitor							x
52		Evaporating temperature setpoint	Control							x
53			Monitor							x
54		Condensing temperature setpoint	Control							x
55			Monitor							x

D: Daikin; ME: Mitsubishi Electric; T: Toshiba Carrier; P: Panasonic.

JH: Johnson Controls-Hitachi Air Conditioning, MH: Mitsubishi Heavy Industries, EM: Emulator

Each company had an item permitting or prohibiting remote control by the occupants. The items that could be permitted or prohibited included the on/off state, operation mode, set-point temperature, airflow volume, and airflow direction, which varied according to the manufacturer. As explained in a later Section 5.1, because we did not develop a detailed model of how these items were manipulated by the occupants, the emulator was designed such that it could only permit or prohibit changing the set-point temperature.

Forced control of the evaporation and condensation temperatures is not currently allowed by either manufacturer. However, many previous studies have reported its energy saving effectiveness (Lian et al. 2005; Miyata et al. 2019; Yun, Lee, and Kim 2016; Zhang et al. 2018). Therefore, we assume that such control will be permitted in the near future, and have added it to the monitoring and control items of this emulator.

Difficulty in measuring the heating/cooling load is a problem with VRF; however, in recent years, attempts have been made to estimate the load from the compressor speed and refrigerant temperature. Therefore, the emulator also included the heating/cooling load of the indoor units as a monitoring item.

Table 6 lists the emulator’s BACnet objects. The emulator must express reasonable changes in comfort and energy when these items are manipulated by a user. Most of the items can be evaluated using the VRF system model already developed by the authors (Togashi and Satoh 2021). For example, this VRF model can be used to evaluate the impact of changing the evaporation and condensation temperatures on energy consumption. In addition, the changes in the energy consumption pattern of the VRF system under different heat loads for each indoor unit can be determined. However, the model needs to be modified for the following two points.

Table 6. BACnet objects related to VRF system in the emulator.

No.	Object name^†1	Type^†2	Description
1	On/Off (C)	BO	This object is used to start or stop the indoor unit.
2	On/Off (M)	BI	This object is used to monitor on/off status of the indoor unit.
3	Operation mode (C)	MO	This object is used to change the operation mode of the indoor unit. 1: cool, 2: heat, 3: fan
4	Operation mode (M)	MI	This object is used to monitor operation mode of the indoor unit.
5	Room temperature setpoint (C)	AV	This object is used to change the setpoint of the indoor unit.
6	Room temperature setpoint (M)	AI	This object is used to monitor the setpoint of the indoor unit.
7	Room temperature (M)	AI	This object is used to monitor room (return) dry-bulb temperature of the indoor unit.
8	Room humidity (M)	AI	This object is used to monitor room (return) relative humidity of the indoor unit.
9	Fan speed or volume (C)	MO	This object is used to change an indoor unit’s fan speed. 1: Low, 2: Middle, 3: High
10	Fan speed or volume (M)	MI	This object is used to monitor the indoor unit’s fan speed.
11	Airflow direction (C)	MO	This object is used to change the indoor unit’s airflow direction. 1: Horizontal, 2: 22.5deg, 3: 45deg, 4: 67.5deg, 5: Vertical
12	Airflow direction (M)	MI	This object is used to monitor the indoor unit’s airflow direction.
13	Fan electricity (M)	AI	This object is used to monitor fan electricity of indoor unit.
14	Heat load of indoor unit (M)	AI	This object is used to monitor heat load of indoor unit.
15	Remote controller Permit / Prohibit room temperature setpoint (C)	BV	This object is used to permit or prohibit the remote controller to set the indoor unit setpoint.
16	Remote controller Permit / Prohibit room temperature setpoint (M)	BI	This object is used to monitor permit or prohibit status of the indoor unit setpoint.
17	Forced evaporating/condensing temperature control (C)	BO	This object is used to change the forced evaporating/condensing control of VRF system.
18	Forced evaporating/condensing temperature control (M)	BI	This object is used to monitor the forced evaporating/condensing control of VRF system.
19	Evaporating temperature setpoint (C)	AV	This object is used to set the evaporating temperature of VRF system.
20	Evaporating temperature setpoint (M)	AI	This object is used to monitor the evaporating temperature of VRF system.
21	Condensing temperature setpoint (C)	AV	This object is used to set the condensing temperature of VRF system.
22	Condensing temperature setpoint (M)	AI	This object is used to monitor the condensing temperature of VRF system.
23	Fan and compressor electricity (M)	AI	This object is used to monitor fan and compressor electricity of outdoor unit.
24	Heat load of VRF system (M)	AI	This object is used to monitor heat load of whole VRF system.

^†1 C: Control, M: Monitor.

^†2 AI: analog input; AO: Analog output; AV: Analog value; BI: Binary input; BO: Binary output; BV: Binary value.

MI: Multistate input, MO: Multistate output.

First, the fan speed and air flow direction (items 9–12) affect the air distribution in the room but are not represented in the general heat load simulator. Thus, here, we explain how to introduce the jet stream model in Section 4.2.

The permission to use remote controllers (items 14 and 15) affects comfort by influencing the psychology of the occupant. Section 5 presents a comprehensive review on thermal comfort and describes the modeling of the impact of occupant psychology on satisfaction.

As controls related to the VRF system, the BACnet objects in Table 7 were prepared for the ventilation system in the emulator.

Table 7. BACnet objects related to the ventilation system in the emulator.

No.	Object name	Type	Description
1	CO2 Level	AI	This object is used to monitor room (return) CO2 level.
2	On/Off	BO	This object is used to start or stop the ventilation unit.
3	Outdoor air bypass enabled	BO	This object is used to bypass outdoor air without heat exchange.
4	Fan speed or volume	MO	This object is used to change a ventilation fan speed. 1: Low, 2: Middle, 3: High

Time synchronization method

If there are multiple controllers, each controller has its own clock. In a real building, both move in real time; therefore, if they are periodically synchronized, there will not be a large discrepancy in time. However, because the emulator can change the acceleration, it sometimes moves much faster than in real time, and sometimes moves slowly in real time. Therefore, if the controller wants to control the equipment on a certain schedule, then it is necessary to synchronize the time with that of the emulator, which changes its speed.

BACnet provides commands for time synchronization (ASHRAE 2020). A device that manages the date and time (called the master) can notify preregistered devices (called slaves) of the current date and time at any given time. Therefore, one solution is to send a time synchronization command from the emulator (master) to the controller (slave); however, this does not actually work. Even if time can be synchronized at a certain moment, if the acceleration speed is not given simultaneously, the time discrepancy widens from the next moment.

This problem was solved using BACnet’s event-notification feature, which notifies users when the value of a variable changes as an event. This is known as the change in value (COV). Each controller registers itself with the emulator to receive COV event notifications when acceleration changes.

When the controller receives a COV event from the emulator, it receives three values related to acceleration: the acceleration rate (A_cc [-]), real-world date and time when the acceleration was changed (DT_real,bs), and date and time inside the emulator at that time (DT_emu,bs). Time cannot be accurately synchronized by acquiring only the acceleration date and the time in the emulator. This is because owing to the time lag required for communication, at the moment the controller acquires the accelerated time, the emulator reaches a virtual time further ahead.

If the function that outputs the difference between two dates and times in seconds is expressed as f_ds (DT₁, DT₂) and the function that adds a specified second to a date and time is expressed as f_add (DT, sec), the accelerated date and time in the emulator (DT_emu,acc) can be calculated using the following equation: (1) $D{T}_{\mathrm{emu},\mathrm{acc}}=f_{\mathrm{add}}\{D{T}_{\mathrm{emu},\mathrm{bs}},A_{\mathrm{cc}}{f}_{\mathrm{ds}}(D{T}_{\mathrm{real}},D{T}_{\mathrm{real},\mathrm{bs}})\}$(1) Here, DT_real represents the date and time of the day in the real world. If this calculation is performed on the emulator and controller sides, the time can be accurately synchronized as long as there are no deviations in the respective clocks.

Building thermal model

Model of heat load calculations

The simulation model to calculate thermal load of building was developed using the program library described by Togashi. The accuracy of the thermal load calculation of this library was verified using BESTEST (Judkoff and Neymark 1995; Togashi and Tanabe 2009). Mass points were set for each of the zones shown in Figure 2 and each layer of the walls. The heat flow between these mass points were expressed via differential equations and solved using the difference method. The solar radiation through a window was calculated based on the incident angle characteristics of the glass.

Each zone was not thermally independent and reflected the heat exchange through air movement as well as wall-to-wall radiant heat transfer. Therefore, the VRF indoor units were thermally affected by zones other than the zone wherein they were placed. The amount of air circulated between zones was assumed to be fixed at 56 CMH for each area where the zones contacted each other. However, in reality, this is influenced by the direction and volume of air blown from the indoor unit. This is due to the use of a simplified airflow model, as explained in the following section, which constitutes the limitation of this emulator’s ability to replicate reality.

Calculation of jet flow from indoor units

VRF controllers can change the air direction by controlling the direction of the air outlet blades. If the airflow is blown horizontally during heating, the vertical temperature distribution increases, resulting in high suction temperatures in the indoor unit and cold temperatures in the occupied area. However, if the airflow is blown directly downwards, it is likely to directly reach the occupants, and there is a risk of a draft. Therefore, the direction of airflow has a significant impact on comfort.

However, general heat-load calculation simulators do not calculate the detailed airflow in a room, making it difficult to represent the above phenomena. CFD should be used to predict airflow with high accuracy; however, its speed is too slow to be used as an emulator. Therefore, we decided to use a jet stream model, which has a relatively small computational load, to predict the airflow from indoor units.

Modeling of zone and blowing jet stream

Figure 3 shows the model of the zone and blowing jet stream.

Figure 3. Model of zone and blowing jet stream.

The room was divided into two parts at FL + 1,700 mm (standing face height), and the temperatures were calculated separately. A mass point is set at the center of both spaces, and the temperatures are expressed as t_U [°C] and t_L [°C], respectively. The vertical temperature gradient dt/dy was calculated from these temperatures.

It is assumed that the indoor unit blows air from the ceiling surface; however, not all of the heat is given to the upper zone, but is divided between the upper and lower zones depending on the jet flow. In Figure 3, h_U [W] represents the heat supplied to the upper zone and h_L [W] represents the heat supplied to the lower zone.

If the jet reaches FL + 1,000 mm (neck height of the seating position), the wind speed at that point (u_m10 [m/s]) was used to calculate the dissatisfaction rate due to the draft. This is discussed in Section 5.2.

The shape of the supply air jet from the indoor unit is affected by the blowing air velocity u_m(0) [m/s], blowing air temperature t_m(0) [°C], blowing air angle θ(0) [degree], and temperature gradient dt/dy.

Notably, when the temperature gradient dt/dy is negative, natural convection prevails and mixes air at the top and bottom; therefore, a large value should be assigned to the ventilation rate between the upper and lower zones.

Jet flow calculation method

To predict the vertical temperature distribution, the path of air blown from an indoor unit must be calculated as it passes through space under the influence of buoyancy owing to the temperature difference. The model proposed by Togari, Arai, and Miura (1993) can calculate such vertical temperature distribution, and Togari reported the results of accuracy verification in large space.

Togari’s model theoretically follows the model proposed by Kubota (1987). Such jet flow paths are particularly important in large spaces where the vertical temperature distribution is large. However. Kubota’s model is a general physical model that is not limited to large spaces; thus, it can be adopted in this study.

Kubota’s model is based on two main assumptions: (1) the velocity and temperature distributions in the cross section perpendicular to the axis of the jet flow direction are similar, and (2) the jet width increases proportionally with the length of the flow path. Kubota compared the predictions of this model with experimental results presented by Cederwall (1967), Frankel and Cumming (1965), and Fan (1967) to confirm its accuracy.

However, Kubota’s model cannot reproduce the interaction of air outlets from multiple indoor units, because this model predicts only one air outlet in space and the way it proceeds in a two-dimensional plane. A three-dimensional CFD analysis would be required to achieve this; however, at present, the calculations would not be completed at a speed suitable for the emulator’s purposes.

Example of jet stream calculation

Using the model described above, the jet stream during heating was calculated for the room shown in Figure 3. Kubota’s model requires the use of the parameters: slot width H_o [m] and throw constant K_p [-]. Assuming a four-way blowing indoor unit, reading the manufacturer’s technical data, the slot width H_o was fixed at 0.5 m and the throw constant K_p was fixed at 4.0.

The jet path calculated from the outlet is shown in Figure 4.

Figure 4. Jet stream model calculation example (heating).

The reference values for air temperature around the air outlet, supply air temperature, blowing angle, blowing air velocity, and vertical temperature gradient were 25°C, 40°C, 45°, 1.5 m/s, and 1.0 K/m, respectively. Each variable was varied within an appropriate range to check sensitivity. The higher the temperature of the airflow, the closer the angle is to the horizontal, and the lower the air velocity, the more likely it is that the air supply will remain in the upper zone. The smaller the vertical temperature gradient, the warmer the air remains in the upper space, but its effect is relatively small compared to the other factors.

The numbers in parentheses in the figure are the heat distribution ratios to the lower space, R_ht [-], calculated using the following equation with the heat supplied to the upper zone, h_H, and the heat supplied to the lower zone, h_L. (2) $R_{\mathrm{ht}}=\frac{h_L}{h_H+h_L}$(2) As shown in Figure 5, for the heat load calculation, independent ducts were considered extending from the indoor unit to each of the upper and lower zones, and the supply airflow was calculated using R_ht in the above formula. It is assumed that air equal to the airflow supplied to the lower zone flows from the lower zone to the upper zone, and the indoor unit draws air from the upper zone.

Figure 5. Assumption of air flow.

The cooling results are shown in Figure 6. During cooling, the airflow path is not as complicated as that during heating because the heavy-density cold air falls straight.

Figure 6. Jet stream model calculation example (cooling).

Calculation method of heat flow not due to jet stream

Even if the jet does not reach the lower zone, the heat transfer will not be zero. Togari et al. (1993) introduced a temperature diffusion coefficient α_t [m²/h] to describe the heat flow when the effect of the jet is small and estimated its value as 3–4 m²/h. The heat transfer coefficient between layers C_a [W/(m²·K)] can be expressed as follows. (3) $C_a=\frac{{\alpha }_t{c}_{\mathrm{pa}}{\gamma }_a}{3600L_s}$(3) where L_s [m] is the thickness of the layer to be calculated, c_pa [J/(kg·K)] is the isobaric specific heat of air, and γ [kg/m³] is density. If we substitute L_s = 1.35 m as the average layer thickness in this model, c_pa = 1000 J/(kg·K), and γ = 1.2 kg/m³, then C_a will be approximately 0.7–1.0 W/(m²·K). In this model, 1.0 W/(m²·K) is used.

Even if the jet stream from the indoor unit does not reach the lower zone, there is generally some air circulation owing to the heat generated by humans or machines in the lower zone, and this air circulation creates an upward air flow. Therefore, we assumed an air circulation equivalent to 0.1 time/h per volume of space, which is less than 1% of the air supply volume of the indoor unit.

Occupant model

This Section describes a method for calculating the dissatisfaction rate of occupants. Although occupants are dissatisfied by many factors, this study was not aimed at constructing a comprehensive model that reflected the impact of all of these factors. The impact of factors that exhibit a relationship with energy consumption must be properly modeled to express the tradeoff between energy performance and comfort.

For this reason, the dissatisfaction rate was estimated for the following four factors: (1) thermal sensations, (2) drafts caused by unwanted airflow, (3) vertical temperature distribution, and (4) poor cleanliness due to lack of ventilation.

Dissatisfied rate due to thermal preference

After predicting the dissatisfaction rate due to thermal preference for each office worker, the dissatisfaction rate was corrected according to whether the remote controller was allowed for the occupants.

Prediction of thermal sensation for each occupant

Classical thermal indices include PMV (Fanger 1970) and SET* (Gagge, Fobelets, and Berglund 1986), but these indices are predictors of the average human thermal sensation. A major challenge in controlling the indoor environment in actual buildings is that not all office workers do not necessarily have the same thermal sensation as the average of the majority, and each individual is different. Some offices have only hot-sensitive workers, whereas others have both hot- and cold-sensitive workers who stay together in the same space. The question is what type of thermal environment should be provided to such a group of workers with different thermal preferences. Therefore, the emulator must also be able to represent the thermal preferences of each office worker. Therefore, this study introduces the thermal preference model proposed by Langevin, Wen, and Gurian (2013).

Langevin’s model is a comprehensive model that can calculate thermal preferences, acceptability, and sensation. Langevin, Wen, and Gurian (2015) also reported a method to simulate the interaction between occupant behaviour and the building environment.

In this model, each office worker has upper and lower limits of thermal sensation for winter (heating) and summer (cooling). These are expressed as s_{th,low,winter}, s_{th,high,winter}, s_{th,low,summer}, and s_{th,high,summer}, respectively, and shall take seven discrete values ranging from -3 (cold) to +3 (hot). The thermal sensation of an office worker is a random variable and its distribution depends on the PMV value (pmv_f [-]), as expressed in the following equation: (4) $P(\mathrm{ASHVote})=f(\mathrm{pm}{v}_f)$(4) The probability that an office worker will not accept a thermal environment owing to cold temperatures in winter and summer is expressed by the following equation: (5) $\begin{array}{r}P(\mathrm{ASHVote}<s_{\mathrm{th},\mathrm{low},\mathrm{winter}}),P(\mathrm{ASHVote}<s_{\mathrm{th},\mathrm{low},\mathrm{summer}})\end{array}$(5) Similarly, the probability that a worker will not accept a thermal environment owing to the heat in each season is expressed by the following equation: (6) $\begin{array}{r}P(\mathrm{ASHVote}>s_{\mathrm{th},\mathrm{high},\mathrm{winter}}),P(\mathrm{ASHVote}>s_{\mathrm{th},\mathrm{high},\mathrm{summer}})\end{array}$(6) Thus, the probability that the n_oc^th occupant is dissatisfied owing to thermal sensation (DR_noc,th [-]) can be calculated seasonally as follows: (7) $D{R}_{n_{\mathrm{oc}},\mathrm{th}}=\{\begin{array}{ll}P(\mathrm{ASHVote}<s_{\mathrm{th},\mathrm{low},\mathrm{winter}})& (\mathrm{winter})\\ +P(\mathrm{ASHVote}>s_{\mathrm{th},\mathrm{high},\mathrm{winter}})& \\ P(\mathrm{ASHVote}<s_{\mathrm{th},\mathrm{low},\mathrm{summer}})& (\mathrm{summer})\\ +P(\mathrm{ASHVote}>s_{\mathrm{th},\mathrm{high},\mathrm{summer}})\end{array}$(7) When there is N_oc number of occupants in a given zone, the average dissatisfaction rate due to thermal sensation (DR_th [-]) is calculated as follows: (8) $D{R}_{\mathrm{th}}=\frac{1}{N_{\mathrm{oc}}}\sum _{n_{\mathrm{oc}}}^{N_{\mathrm{oc}}}D{R}_{n_{\mathrm{oc}},\mathrm{th}}$(8) Langevin, Wen, and Gurian (2015) also proposed a model for how office workers adjust their clothing, which is adopted in this study. That is, the office worker adjusts the amount of clothing (clo_mrn [clo]) on the morning of the n^th day according to the n − 1^th day’s thermal environment experience using the following equation: (9) $\begin{array}{rl}\log (\mathrm{cl}{o}_{\mathrm{mrn}}(n))& =-0.91-0.01t_{\mathrm{out},6\mathrm{am}}\\ & +0.14W_{\mathrm{pref}}+0.71\mathrm{cl}{o}_{\mathrm{mrn}}(n-1)\end{array}$(9) where t_out,6am [°C] is the outdoor air temperature at 6 am, and W_pref [-] is a dummy variable that takes the value of 1 if the thermal neutrality of the occupant is on the warm side. In addition, the occupant was given the opportunity to adjust the amount of clothing by rolling up the sleeves (0.08 clo) and adding or removing the sweater (0.3 clo) at every time step (15 min). If the occupant is not satisfied with the clothing adjustment, he/she tries to adjust the temperature using the HVAC remote controller.

The information related to occupants that can be monitored by BACnet communication includes the presence/absence of each occupant, the seven levels of thermal sensation calculated using the method described above, and the amount of clothing worn. Although this information is not currently available for real buildings, it is expected to be possible in the future through computerization.

Correction by remote control permission

A major characteristic of VRF systems compared to general central heat source systems is that the occupant is usually allowed to manipulate the setpoint temperature and other settings. The results of previous studies, such as the following, suggest that such a feature is likely to have an impact on the satisfaction of the occupants.

The same physical stimulus does not always produce the same sensation in a person; rather, the sensation varies depending on the person’s state. This phenomenon was defined by Cabanac, who is psychophysiologist, with the term ‘Alliesthesia’ (Cabanac 1999). Subsequently, Parkinson and de Dear (2015) applied this term to the hypothesis of ‘thermal alliesthesia’ and summarized it as follows: That is, ‘Any peripheral thermal stimulus that offsets or counters a thermoregulatory load-error will be pleasantly perceived’.

In conventional models, such as PMV and SET*, a person is considered a static machine that outputs the thermal sensation only according to the instantaneous values of the six thermal factors. However, according to the hypothesis of thermal alliesthesia, human sensation is influenced by the thermal history of a time series, making people proactive entities who try to act in future thermal environments.

In fact, several studies have reported that manipulating a person’s thermal environment changes the person’s perception of the environment more than the physical changes brought about by the manipulation. For instance, Boerstra et al. (2013) analyzed the HOPE database and found a statistically significant relationship between the degree to which people felt they had control over their thermal environment and the degree to which they felt comfortable in that thermal environment. Kwon et al. (2019), Brager, Paliaga, and de Dear (2004), and Chen et al. reported similar findings. Through laboratory experiments, Schweiker et al. (2012) reported that the ability to operate tabletop fans and open and closed windows increased satisfaction.

An early example of an attempt to model the impact of such controllability on thermal comfort was Paciuk (1989), who defined the key concepts of ‘Available control’, ‘Exercised control’, and ‘Perceived control’. Hellwig (2015) proposed a comprehensive model for the relationship between personal control and thermal sensation. However, these conceptual models do not have specific formulae or parameters. Ackerly, Brager, and Arens (2012) proposed standard measurement and questionnaire methods for perceived control, and Hawighorst, Schweiker, and Wagner (2016). Thus, data are collected; however, it is expected that it will take some time before a generic model is developed.

Therefore, the model to be implemented in an emulator cannot be comprehensive, and must be simplified. Currently, the concept of ‘adaptive increments’ (Baker and Standeven 1996) would be easier to implement. People expand their range of comfort by proactively controlling and adapting to their environment, and adaptive increments express this magnitude in terms of temperature units.

For example, Yun (2018) studied seven university facilities in South Korea and reported that when thermostat control was allowed for office workers, the comfort temperature during cooling rose from 24.3 to 25.2°C. Although not temperature control, Haldi and Robinson (2008) studied the values of adaptive increments when window openings, blind openings, and the use of tabletop fans were allowed and reported values of 0.58, 0.46, and 1.39°C, respectively. In a study on multifamily buildings rather than commercial buildings, Luo et al. (2014) reported a 2.6°C, change in thermally neutral temperatures when controllers are allowed during heating.

From the above, it can say that allowing control is roughly equivalent to changing the temperature by about 1 °C, which is about 0.2 in PMV. This is expressed in Eq. (4) to express the effect of allowing personal control on reducing dissatisfaction. Specifically, we first calculated the value of the thermal sensation that each occupant felt thermally neutral using the following equation: (10) $\begin{array}{rl}{s}_{\mathrm{th},\mathrm{neutral},\mathrm{winter}}& =0.5(s_{\mathrm{th},\mathrm{low},\mathrm{winter}}+s_{\mathrm{th},\mathrm{high},\mathrm{winter}})\end{array}$(10) (11) $\begin{array}{rl}{s}_{\mathrm{th},\mathrm{neutral},\mathrm{summer}}& =0.5(s_{\mathrm{th},\mathrm{low},\mathrm{summer}}+s_{\mathrm{th},\mathrm{high},\mathrm{summer}})\end{array}$(11) According to the following equation, we can calculate the corrected PMV value (pmv_cr [-]) using the above thermally neutral sensation and Fanger’s PMV value, and substitute this for pmv_f in Eq. (4). (12) $\mathrm{pm}{v}_{\mathrm{cr}}=\{\begin{array}{ll}\mathrm{pm}{v}_f+min(c_{\mathrm{tstt}},& (\mathrm{pm}{v}_f<s_{\mathrm{th},\mathrm{neutral}})\\ s_{\mathrm{th},\mathrm{neutral}}-\mathrm{pm}{v}_f)& \\ \mathrm{pm}{v}_f-min(c_{\mathrm{tstt}},& (s_{\mathrm{th},\mathrm{neutral}}<\mathrm{pm}{v}_f)\\ \mathrm{pm}{v}_f-s_{\mathrm{th},\mathrm{neutral}})\end{array}$(12) where c_tstt [-] is the PMV correction value, which is 0.2 when the occupants are allowed to control their own environment and 0.0 when environment control is prohibited.

Dissatisfied rate due to draft

According to ISO 7730 (2005), the dissatisfaction rate DR_dr,ISO [-] owing to the draft can be predicted using the following equation: (13) $\begin{array}{rl}D{R}_{\mathrm{dr},\mathrm{ISO}}& =(34-t_{\mathrm{lcl}})(\overline{u_{\mathrm{lcl}}}-0.05{)}^{0.62}(0.0037\overline{u_{\mathrm{lcl}}}{T}_{\mathrm{urb}}\\ & +0.0314)\end{array}$(13) where t_lcl [°C] is the local air temperature; in ISO7730, the applicable range of the equation is 20–26°C, but in this model, the range is up to 34°C in order to reflect differences in blow-off temperature in the results. $\overline{u_{\mathrm{lcl}}}$ [m/s] is the local mean air velocity, $\overline{u_{\mathrm{lcl}}}$ = 0.05 when $\overline{u_{\mathrm{lcl}}}$ < 0.05. T_urb [-] is the local turbulence intensity, which can be calculated using the standard deviation of the local airflow, δ_lcl [m/s], as follows. (14) $T_{\mathrm{urb}}=\frac{{\delta }_{\mathrm{lcl}}}{\overline{u_{\mathrm{lcl}}}}$(14) T_urb [-] takes values from 10% to 60%; however, when unknown, 40% is allowed, which was fixed at 40% in this model.

While CFD, which solves for the airflow in the entire 3D space of a room, can provide information on the velocity distribution, this model solves for the jet stream flow path in a 2D plane; therefore, the calculation results cannot be directly applied to Eq. (13). Therefore, as shown in Figure 7, the space was divided into concentric circles centered on the jet stream, and the following procedure was used to predict the number of dissatisfied occupants in the room.

Figure 7. Dividing space for draft calculations.

In the Kubota model described earlier, the velocity u [m/s] and temperature t [°C] at a distance r [m] perpendicular to the axis from the center of the jet are expressed as follows: (15) $\begin{array}{rl}u& =u_m\exp \left\{-\frac{\pi }{2}{\left(\frac{K_{p}r}{s}\right)}^2\right\}\end{array}$(15) (16) $\begin{array}{rl}t& =t_e-(t_e-t_m)\exp \left\{-\lambda \frac{\pi }{2}{\left(\frac{K_{p}r}{s}\right)}^2\right\}\end{array}$(16) where s [m] is the distance along the jet from the air outlet, u_m [m/s] and t_m [°C] are the velocity and temperature of the jet center, respectively, and t_e [°C] is the ambient temperature.

Solving Eq. (15) for r, calculate the distance r_dr [m] at which a draft can occur as follows: (17) $r_{\mathrm{dr}}=\frac{s}{K_p}\sqrt{\frac{2}{\pi }\ln \left(\frac{u_m}{u_{\mathrm{dr},\min }}\right)}$(17) where u_dr,min [m/s] is the minimum velocity at which a draft can occur, which is 0.05 m/s on the definition of Eq. (13).

This r_dr is divided equally into N_r and the distance from the center of the jet is expressed as r_dr,0= 0, r_dr,1, r_dr,2 … r_dr,Nr. That is, a circle of radius r_dr is divided into N_r annuli of width r_dr/N_r.

The central diameter r_an,n [m] of each annulus was calculated using the following equation: (18) $r_{\mathrm{an},n}=0.5(r_{\mathrm{dr},n}+r_{\mathrm{dr},n+1})$(18) The values obtained by substituting r_an,n into Eq. (15) and Eq. (16) are regarded as the representative wind speeds, u_an,n and representative temperatures, t_an,n in each annulus area. These are then substituted into Eq. (13) to estimate the dissatisfied rate DR_an,n [-] for each annulus area. The overall dissatisfaction rate by the jet stream DR_jet [-] is estimated by weighted averaging of these dissatisfaction rates DR_an,n by the area of the annulus as follows: (19) $D{R}_{\mathrm{jet}}=\sum _{n=0}^{N_r-1}{w}_{\mathrm{an},n}D{R}_{\mathrm{an},n}$(19) The weighting factor w_an,n [-] is calculated by follows. (20) $w_{\mathrm{an},n}=\frac{r_{\mathrm{an},n+1}^2-r_{\mathrm{an},n}^2}{r_{\mathrm{dr}}^2}$(20) DR_jet is the probability of dissatisfaction when exposed to a jet stream, whereas no dissatisfaction occurs when the occupants are seated in a seat that is not exposed to the jet stream. The three-dimensional relationship between the jet stream and the occupant’s seat is necessary to determine whether the occupant is exposed to the jet; however, this is not represented in the emulator. However, on average, if there are more occupants in a zone, the likelihood of being hit by a jet stream is higher, and if there are fewer occupants, the likelihood is lower. Thus, using the occupancy rate R_oc,c [-] (values: 0.0–1.0), the dissatisfied rate due to draft for the entire zone (DR_drft [-]) is predicted by the following equation. (21) $D{R}_{\mathrm{drft}}=R_{\mathrm{oc},c}D{R}_{\mathrm{jet}}$(21) Such an assumption implies that it is effective to air condition the room strongly in advance in the morning or during lunch breaks when there are relatively few occupants, and to blow the air as gently as possible during hours when there are many occupants, which would also be a viable practice in reality.

Notably, drafts are mainly a problem when they are cold, whereas a moderate airflow is pleasant when they are hot. Therefore, the occupancy rate R_oc,c was calculated only for occupants whose thermal sensations were on the hot side (slightly warm, warm, hot).

Figure 8 shows the calculation results for the variation in the draft dissatisfaction rate DR_jet with the combination of the velocity and temperature of the airflow. The temperature t_m10 and velocity u_m10 at the center of the jet at a height of 1.0 m and the ratio of heat supplied to the lower zone, R_ht, are also shown.

Figure 8. Dissatisfied rate due to draft (DR_jet) relative to supply air velocity and temperature.

The velocity of the jet center was high at high blowing velocities, whereas it was low at high blowing temperatures owing to velocity reduction due to buoyancy. By contrast, the temperature at the center of the jet stream was high at high blowing temperatures. Moreover, the temperature is slightly high at high blowing velocities because the time required to reach a height of 1.0 m is less, and less heat is dissipated from the jet to the surrounding air.

According to the above trend, the percentage of dissatisfied occupants increased with a lower blowing temperature and higher blowing velocity, that is toward the lower right of the figure. Conversely, in the area above the dotted line in the figure, the temperature is high and the velocity is low; therefore, the airflow does not reach the occupant. In this case, no draft can occur, and the dissatisfaction ratio is zero.

Therefore, to reduce the risk of drafts, indoor units should be operated in the upper-left region of the figure. However, this operation is problematic in terms of temperature control in the zone, because it is difficult for hot air to reach the lower space. Looking at the diagram of the heat supply ratio to the lower zone (R_ht), the upper left panel of the diagram shows smaller values. In other words, a heat pool occurs at the top, resulting in poor performance in terms of vertical temperature distribution.

In this calculation example, there is a region where the ratio of heat supplied to the lower zone (R_ht) is not zero, whereas the dissatisfaction rate owing to the draft is kept at zero. If the indoor unit is operated in this region, neither vertical temperature differences nor drafts occur. However, in reality, it is difficult to continue operation in such regions because the room temperature changes continuously. As the thermal capacity of a building allows heat to be stored, it may be effective to use air conditioning more intensively when there are fewer occupants, as mentioned earlier.

In terms of the energy efficiency of the machine, changes in the blowing temperature affect the efficiency of the refrigeration cycle through the condensation or evaporation temperature, and the blowing air speed affects the heat exchange efficiency between the refrigerant and air in the indoor unit. To reduce draft risk, the blowing angle can be adjusted in addition to changing the blowing air speed. In terms of heat exchange capacity, the former, which can maintain airflow, is more advantageous and increases the energy efficiency of the machine.

In other words, we must solve the difficult problem of optimizing the operation from three perspectives: draft risk, vertical temperature distribution, and machine energy efficiency.

Dissatisfied rate due to vertical temperature difference

In accordance with ISO 7730 (2005), the dissatisfaction rate DR_vtmp [-] due to the vertical temperature difference is predicted using the following equation: (22) $D{R}_{\mathrm{vtmp}}=\frac{1.0}{1.0+\exp (5.76-0.856\mathrm{\Delta }{t}_{\mathrm{vtmp}})}$(22) where Δt_vtmp [°C] is the difference in temperature at foot and head height (0.1 m and 1.1 m), obtained by linear interpolation using the temperatures at the upper and lower zones calculated by the model described in section ラー! 参照元が見つかりません。.

Dissatisfied rate due to lack of ventilation

Generally, the cleanliness of indoor air is estimated based on the concentration of carbon dioxide. For example, Japanese law sets an upper limit of 1,000 ppm for many building applications. Since carbon dioxide concentration is an indirect indicator, people are not necessarily able to recognize it even if the concentration exceeds 1,000 ppm. However, in this model, the following equation is used to generate dissatisfaction: (23) $D{R}_{\mathrm{co}2}=\frac{1.0}{\exp \left(\frac{1050-L_{\mathrm{co}2}}{10.8}\right)}$(23) This was based on the assumption that 1% dissatisfaction occurred at 1000 ppm and 99% at 1100 ppm, as shown in Figure 9. The purpose of this assumption is not to correctly predict dissatisfaction rates but to penalize operations that reduce ventilation too much to reduce energy.

Figure 9. Dissatisfied rate according to CO2 level.

Integration of the four dissatisfied rates

If each dissatisfaction calculated above is considered to occur independently, the total dissatisfaction rate DR_ttl [-] can be calculated as follows: (24) $\begin{array}{rl}D{R}_{\mathrm{ttl}}& =1-(1-D{R}_{\mathrm{th}})(1-D{R}_{\mathrm{drft}})\\ & \times (1-D{R}_{\mathrm{vtmp}})(1-D{R}_{\mathrm{co}2})\end{array}$(24) Because DR_zn varies with time and zone, it is weighted by the number of occupants in that zone at that time step to obtain an average value, which is the average total dissatisfaction rate, DR [-], for the entire calculation period. In other words, if the number of zones is N_zn [-], and the dissatisfaction rate and number of occupants in each zone at time step i are DR_ttl,zn,i [-] and N_oc,zn,i [persons], respectively, the average dissatisfaction rate DR can be calculated as follows: (25) $\mathrm{DR}=\frac{{\sum }_i{\sum }_{\mathrm{zn}}^{N_{\mathrm{zn}}}(N_{\mathrm{zn},i}D{R}_{\mathrm{ttl},\mathrm{zn},i})}{{\sum }_i{\sum }_{\mathrm{zn}}^{N_{\mathrm{zn}}}{N}_{\mathrm{zn},i}}$(25)

Consideration of championship rules

This section presents the results of several rule changes based on the reflections of the first WCCBO.

Total score calculation method

In the first WCCBO, the total score (E2R2 [-]) was calculated using the following formula: (26) $E2R2=\{\begin{array}{ll}\mathrm{ERR}\cdot \mathrm{DRR}& (0<\mathrm{ERR}\cap 0<\mathrm{DRR})\\ 0& (\mathrm{other})\end{array}$(26) where ERR [-] (Energy Reduction Ratio) and DRR [-] (Dissatisfaction Reduction Ratio) are calculated by the following equations, in which E [MJ] is the primary energy consumption, D [persons] is the total number of dissatisfied occupants, the subscript r indicates reference operation, and the subscript opt indicates optimized operation. (27) $\begin{array}{rl}\mathrm{ERR}& =\frac{E_r-E_{\mathrm{opt}}}{E_r}\end{array}$(27) (28) $\begin{array}{rl}\mathrm{DRR}& =\frac{D_r-D_{\mathrm{opt}}}{D_r}\end{array}$(28) Equation (26) is established to balance energy conservation and comfort during the operation. The equation clearly shows that the score for an operation that performs poorly in either of these aspects is zero.

However, there was no definitive method to determine reference energy consumption E_r and reference dissatisfaction rate D_r. In the first WCCBO, the E_r and D_r values were determined by testing various operating parameters in the committee prior to the championship and assuming typical operations. In other words, typical operations and reference values are not scientific, but vary depending on the committee’s considerations.

In fact, the first WCCBO used a smaller energy reduction baseline, making it more difficult to improve energy consumption than comfort. The distribution of the competitors’ scores under these conditions is shown in Figure 10 (left). The horizontal axis represents the energy reduction ratio, and the vertical axis represents the unsatisfactory reduction ratio. The horizontal distribution was narrower owing to the difficulty in energy reduction. However, the distribution in the vertical direction was wider, because there was much room for improvement in comfort.

Figure 10. Distribution of scores when it is more difficult to improve energy performance than comfort.

If such a distribution is the case, for example, on the right side of Figure 10, a move from point (a) to point (b) owing to a changing operation should be rated higher than a move from point (a) to point (c). Although the moving distance was the same, the former reached an unexplored territory and clearly solved a difficult problem. In other words, those who could accomplish a relatively difficult operation were rated higher.

However, as Eq. (26) expresses a ratio to a reference value, a reduction of one unit of energy and one unit of dissatisfaction would be evaluated as the same score.

In light of these issues, the second WCCBO considered a method for calculating scores using the distribution of competitors’ operational results. This is a dynamic evaluation method, because the distribution changes each time a new operational result is added.

The simplest approach would be to simply transform the datasets of energy consumption and dissatisfied number into a standard normal distribution, but this is not appropriate. For the 339 operation results of the first WCCBO, the frequency distributions of energy consumption and dissatisfied number are shown in Figure 11. When tested for normality, the original data were not normally distributed.

Figure 11. Frequency distribution of energy consumption and dissatisfied number.

A major reason for normality that cannot be affirmed is the assumption due to the asymmetry in the distribution. In general, although there are limits to reducing energy or increasing comfort, it is easy to increase energy or decrease comfort, resulting in large outliers. For example, to reduce the energy consumption by 50% via the appropriate operation of an HVAC is clearly difficult than increasing it by 50% by operating it poorly.

Figure 12 shows a scatter plot of the number of dissatisfiers and energy consumption. Outliers can be eliminated by removing the data using lower 3% or 6% thresholds.

Figure 12. Scatter plot of the number of dissatisfiers and energy consumption.

In addition, as this was a simulation, identical input conditions produced the same calculation results. Because allowing these duplicates would distort the distribution, data that were exactly the same were excluded.

The processed data were Box–Cox transformed to further approximate a normal distribution and then transformed to a standard normal distribution with a mean of 0 and standard deviation of 1. The frequency distribution of the transformed data is shown in Figure 13. The results of calculating p-values for the normality tests are listed in Table 8. Although imperfect, normality was confirmed by several tests.

Figure 13. Result of transformation to standard normal distribution.

Table 8. Normality test results (p-value).

Test method	Before Box-Cox transformation		After Box-Cox transformation
	Dissatisfied number	Energy consumption	Dissatisfied number	Energy consumption
Shapiro-Francia	0.010	0.000	0.042	0.000
D’Agostino-Pearson	0.166	0.000	0.547	0.000
Jarque-Bera	0.190	0.001	0.627	0.033
Cramer-von Mises	0.003	0.000	0.018	0.100
Anderson-Darling	0.003	0.000	0.018	0.036

† under lined: p = 0.05.

Thus far, we have considered using the total dissatisfied number D to use the data from the first WCCBO, but there is no problem with using the average total dissatisfied rate DR defined in Eq. (25) instead. After transforming the energy consumption E and average total dissatisfaction rate DR into a standard normal distribution, the values expressed in percentiles were expressed as P_E [-] and P_DR [-], respectively. For example, P_E = 90% implies the top 10% score. The ECI (Energy efficiency and Comfort Index (ECI)), an integrated score of energy reduction and comfort, was defined using the following equation: (29) $\mathrm{ECI}=\sqrt{P_E{P}_{\mathrm{DR}}}$(29)

Network structure for score registration

In the first WCCBO, the emulators were placed on a server provided by the organizers and championship participants controlled the server-side emulators online using their local computers. The objective was to familiarize the participants with such remote controls, as it can be assumed that most buildings will be optimized from the cloud in the near future.

According to the difficulty level of remote optimization, the championship was divided into the first and second halves. In the first half, many people participated because the first half was a simple method of sending a batch of schedules entered on a spreadsheet software sheet to the emulator. The second half of the competition, on the other hand, was a difficult method that involved applying VPN technology to connect to an emulator on a server and change control values in real time. Because all settings could be changed, the participants were able to consider more advanced controls, aiming to increase the difficulty level as the competition progressed.

However, according to the post-championship hearings, few teams participated in real-time control in the second half of the championship, and only one team wrote its own optimization program. Unfortunately, the team could not improve their results using real-time control.

In light of the above, it can be concluded that remote optimization control is still too difficult for many in our field, and that it is too early to meet this requirement. Therefore, in the second WCCBO, the emulator was placed on a local PC. The structure of the network used to register the scores is illustrated in Figure 14.

Figure 14. Network structure to register score.

The red frame represents the range in which the organizers of the championship prepare and provide the participants. The emulator ran on the participant’s personal computer (PC). BACnet devices exist inside the emulator, and the BACnet controller inside the local PC or the BACnet controller connected to the network to which the local PC belongs can operate these devices from outside the emulator.

A secret key was embedded inside the emulator at the source code level, and the calculation results were encrypted using the ChaCha20-Ply305 method (Nir and Langley 2015). When a championship participant uploaded the result file to the server, the data were saved after decryption with the server-side key, and the scores were published on the web server.

Although this method can prevent data tampering during communication, it cannot completely prevent data tampering by championship participants themselves who use emulators. Because the emulator is placed on the participant’s local PC, even if the secret key is directly embedded in the source code and compiled, in principle, a reverse compilation cannot be prevented and the secret key will be exposed. However, this degree of robustness would be sufficient as there would be little incentive to falsify data to obtain an unexplained high score if the participants with the best scores were obligated to explain what controls they had attempted.

How to express uncertainty

In a typical building energy simulation, the room usage and outdoor air conditions are entered as fixed boundary conditions. However, these conditions are stochastic and cannot be predicted in advance, making it difficult to operate an HVAC system. To express this point, the first WCCBO introduced a stochastic model (Togashi 2018) for the behaviour of occupants and weather conditions during both the first and second halves of the championship. However, during the first half of the championship, the seed for generating random numbers was fixed so that participants could know inherently stochastic information in advance. However, this does not truly assess the skill of HVAC system operation.

When we obtain data by experiencing a stochastic phenomenon, it is not one specific instance of the data that we should be learning, but the parameters of the stochastic model behind the phenomenon, because if we overadapt to data that are only one instance, we will be powerless against unknown data. In machine learning, this is known as the problem of overfitting. In recent years, machine learning has been increasingly applied in the field of building equipment, and it is important to know how to represent uncertainty in an emulator to use it as a testbed for these studies.

The emulator for the second WCCBO provided two random number seeds to represent the uncertainty associated with the occupants.

Based on the random numbers produced by the first random seed, the thermal preferences (hot- or cold-sensitive) of the occupants in each zone and the behaviour patterns they tend to adopt are initialized. These are the parameters of the stochastic model and the information that must be learned.

According to the random numbers produced by the second random seed, the office worker expresses a thermal sensation or enters or leaves the office. These are only examples of non-reproducible realizations, but they can be used as clues to the parameters of the background stochastic model.

Depending on the values of these random seeds, the championship was divided into two periods: preparation and evaluation. Figure 15 shows the use of random seeds.

Figure 15. Method for utilizing random seed.

The first random number of seeds was fixed during the preparation and evaluation periods. The parameters generated by this seed must not change because these parameters are the values that are to be estimated. A single calculation yields one week of data; however, provided the value of this seed is not changed, we can generate as much training data as required by repeating the calculation.

However, the second random seed may or may not be fixed during the preparation period. If this value was fixed, the results of the case study would be easier to interpret because the emulator would remove the randomness of the weather conditions and occupant behaviour. However, optimizing the operation based on such a case study may result in overfitting. Thus, to avoid overfitting, the user should not fix the second seed and observe if the optimized operation method performs equally well for other weather conditions and occupant behaviour.

At the beginning of the evaluation period, a second random seed value was specified, and the score when this value was used determined ranking. This eliminates the randomness caused by random seed and represents only the differences in skill of HVAC system operation. If the participants have sufficiently identified the parameters of the stochastic model during the preparation period, they will perform well using the same control methods as during the preparation period. Conversely, if the control insights obtained during the preparation period were merely overfitting, they would not be able to cope with the new random number seed value and would be ranked lower. Therefore, the length of the evaluation period should be sufficiently short so that participants do not have time to overfit their optimization models again with the value of the second random seed disclosed.

Sensitivity analysis

List of simulation cases

Several simulations were performed under different conditions using the emulator developed in the previous sections. Table 9 presents the simulation results.

Table 9. Conditions of simulation cases.

Case	-	Condensing / Evaporating temperature [°C]	Setpoint temperature [°C]	Fan speed^†	Airflow direction [degree]	Remote control permission	Stop VRF in the interior zone
H1	Heating	46.0	22.0	Middle	45.0	False	False
H2		40.0	22.0	Middle	45.0	False	False
H3		46.0	26.0	Middle	45.0	False	False
H4		46.0	22.0	Low	45.0	False	False
H5		46.0	22.0	Middle	5.0	False	False
H6		46.0	22.0	Middle	90.0	False	False
H7		46.0	22.0	Middle	45.0	True	False
H8		46.0	22.0	Middle	45.0	False	True
C1	Cooling	10.0	26.0	Middle	45.0	False	False
C2		15.0	26.0	Middle	45.0	False	False
C3		10.0	22.0	Middle	45.0	False	False
C4		10.0	26.0	Low	45.0	False	False
C5		10.0	26.0	Middle	5.0	False	False
C6		10.0	26.0	Middle	90.0	False	False
C7		10.0	26.0	Middle	45.0	True	False
C8		10.0	26.0	Middle	45.0	False	True

^† Middle and Low is 70% and 30% air flow rate to nominal air flow rate, respectively.

H1–H8 were calculated for winter heating, and C1–C8 for summer cooling. H1 and C1 are the standard reference conditions, and H2–H7 and C2–C7 are calculated by changing the refrigerant temperature, room setpoint temperature, fan speed, airflow direction, and remote control permission, respectively. In addition, H8 and C8 were established as cases in which the VRF in the interior zone was forcibly stopped, and air conditioning was performed only with the VRF in the perimeter zone.

The calculation results are shown in Table 10, 11, and Figure 16. The P_E, P_D, and ECI values were calculated using Eq. (29) by using data from eight cases of heating and cooling.

Figure 16. Simulation results.

Table 10. Simulation results (Heating cases).

Case	H1	H2	H3	H4	H5	H6	H7	H8
Energy consumption [GJ]	5.190	4.791	7.432	4.891	4.727	5.281	7.503	4.299
Averaged dissatisfaction rate [-]	0.249	0.256	0.209	0.246	0.238	0.245	0.188	0.222
PE [-]	0.438	0.674	0.098	0.608	0.718	0.397	0.096	0.956
PDR [-]	0.214	0.068	0.864	0.295	0.520	0.329	0.916	0.769
ECI	0.306	0.214	0.291	0.424	0.611	0.361	0.297	0.858

Table 11. Simulation results (Cooling cases).

Case	C1	C2	C3	C4	C5	C6	C7	C8
Energy consumption [GJ]	8.156	6.830	9.296	8.152	8.154	8.158	9.262	6.823
Averaged dissatisfaction rate [-]	0.222	0.226	0.227	0.215	0.207	0.227	0.213	0.251
PE [-]	0.501	0.906	0.089	0.503	0.502	0.500	0.096	0.907
PDR [-]	0.456	0.346	0.320	0.711	0.950	0.330	0.804	0.051
ECI	0.478	0.560	0.169	0.598	0.690	0.406	0.278	0.214

A case-by-case discussion is reported below. However, in general, the dissatisfied rate tends to be larger in the heating case (Hn). Occupants were in the lower part of the room, but warm air tends to accumulate near the ceiling. However, when high winds are used to deliver warm air to the lower part of the room, the potential for drafts increases. These tradeoffs inherent in heating case resulted in the above tendencies.

Change in refrigerant temperature (Case H2, C2)

Lowering the condensing temperature of the indoor unit during heating or raising the evaporating temperature during cooling reduces energy consumption because the compression ratio of the refrigeration cycle is lowered. Therefore, H2 and C2 consumed less energy than H1 and C1. The relationship between the partial-load rate and electricity for each case in the VRF1 system is shown in Figure 17.

Figure 17. Relation between partial load rate and electricity of VRF1. (Case H1 and H2, C1 and C2).

As the figure clearly shows, even at the same partial load rate, the electricity consumption is lower when the evaporation and condensation temperatures are changed. There was an operation with extremely high electricity at a partial load rate of 40–50% during heating, which was caused by a defrosting operation at dawn when the outdoor air temperature was close to 0°C.

As shown in Figure 17, there are a few hours in which the electricity reaches its maximum value and becomes overloaded. In such cases, even if the condensing or evaporating temperature is changed, the indoor thermal environment does not deteriorate to an extreme degree because the capacity of the indoor unit is sufficient and there is little impact on the dissatisfaction rate.

However, care should be taken during the cooling operations. Figure 18 shows the temperature and humidity during business hours in S1. In case C2, where the evaporating temperature was raised, the room humidity increased because the supply air temperature was higher, and the latent heat rate decreased. Although this has little impact on thermal comfort, extreme temperatures can deteriorate a room’s thermal environment.

Figure 18. Air state of zone S1 (Case C1 and C2).

Change in setpoint temperature (Case H3, C3)

If the room setpoint temperature is raised during heating (H1–H3), then energy consumption increases, comfort improves, and the dissatisfaction rate decreases. Conversely, lowering the room setpoint temperature during cooling (C1→C3) only increases the energy consumption and hardly changes the dissatisfied rate because of the following reasons.

The ratios of the thermal sensations declared by the occupants are shown in Figure 19. It is clear that lowering the room setpoint temperature improves the thermal environment because the majority of the thermal sensations have moved to neutral.

Figure 19. Percentage of total thermal sensation of occupants (Case C1 and C3).

Figure 20 shows the change in the dissatisfaction rate on a representative day. Figure 19 indicates that the dissatisfaction rate due to thermal sensation is lower than that of the reference case C1, whereas the rate due to the draft is higher than that of C1. This result can be attributed to the low temperature setting, which increases the number of occupants with a thermal sensation on the cooler side, making them more sensitive to drafts. These results suggest that simply lowering the temperature is not sufficient to make the occupants comfortable during the summer and that the optimal setpoint temperature at which the total dissatisfaction rate is minimized must be determined by considering the trade-off between thermal sensation and drafts.

Figure 20. Time changes in the dissatisfied rate (Case C1 and C3).

Both case exhibited a large increase in the dissatisfied rate in the evening. This was attributed to the air conditioning being turned off in the evening. However, as only a few occupants remain in the office in the evening owing to overtime work, this effect is small when the dissatisfied occupant rate is calculated on a daily average basis.

Change in fan speed and airflow direction (cases H4, H5, H6, C4, C5, and C6)

Lowering the fan speed (H4) or moving the air closer to the horizontal direction (H5) during heating lowers energy consumption. This was because it became difficult for the blown airflow to reach the occupied zone, and the temperature in the lower zone decreased. Therefore, the heating load of the indoor units decreases and energy consumption decreases. If the airflow direction is vertical (H6), the hot air reaches the lower zone. The energy consumption in this case was almost the same as that in the reference case.

Although the dissatisfaction rates for these cases appear to be nearly equal, the actual breakdowns differ. Figure 21 shows the lower zone temperature at the start of air conditioning. Compared with the reference case (H1), H4 and H5 had less warm air reaching the lower zone, resulting in a slower temperature increase. Therefore, the thermal sensation of the occupants shifts to the cold side, resulting in greater dissatisfaction.

Figure 21. Lower zone temperature at start of air conditioning (H1 vs H4, H5).

Figure 22 shows the dissatisfaction rate owing to the draft for each case relative to reference case H1. Cases H4 and H5 had less dissatisfaction owing to the draft, because the airflow had difficulty reaching the lower zone. Case H6, was almost the same as the reference case.

Figure 22 Comparison of dissatisfied rate (H1 vs H4, H5, H6).

In other words, in Cases H4 and H5, the dissatisfaction rate due to thermal sensation increases, but the dissatisfaction rate due to the draft decreases, and the combined effect of the two is that the dissatisfaction rate is almost the same as that in the reference case.

Lowering the fan speed (C4) or moving the air direction closer to the horizontal direction (C5) during cooling would slightly reduce the dissatisfaction rate compared to blowing vertically (C6). This is because cold air sinks naturally and the possibility of drafting is lower when blowing horizontally and gently.

Change in remote control permission (Case H7, C7)

When remote control was allowed, the dissatisfaction rate decreased for both heating (H7) and cooling (C7). In contrast, the energy consumption increased in both cases.

Figure 23 shows the temperature set points of the VRF1 system for one week. The thin line represents the setpoint of each indoor unit and the thick line represents the average. The setpoints of all indoor units do not necessarily coincide with each other because the thermal preferences and amount of clothing worn differ for each occupant. However, the major trend, including other VRFs, is that the setpoint temperature is raised slightly for heating compared with the initial setpoint (22°C for heating and 26°C for cooling) and conversely decreased slightly for cooling.

Figure 23. Time changes in the setpoint temperature (H7, C7).

Notably, although certain zones are set to 32°C for heating or 16°C for cooling, this does not imply that the occupant really needs that temperature. The occupant will only raise or lower the setpoint according to his/her thermal perception (i.e. hot or cold). For example, if a hot worker who feels comfortable at 20°C is in the 21°C zone, he or she will continue to lower the setpoint, such that the setpoint eventually reaches the lower limit of 16°C. As explained earlier, as the zones are thermally connected to each other to a certain extent by air movement, even if the setpoint of a particular zone is only lowered, the setpoint is not easily realized because relatively warm air flows in from the neighbouring zones. Consequently, the setpoints of zones with hot or cold occupants tend to exhibit extreme values.

Figure 24 shows the average clothing worn by the office workers. For heating, the weekly amount of clothing was 0.946 in the reference case (H1) when remote control was prohibited and dropped to 0.838 when remote control was permitted. In other words, while occupants set the room temperature to a higher level with the remote control, they lowered their clothing. Conversely, for cooling, the temperature is lowered by remote control, while on average, more clothes are worn than in the reference case.

Figure 24. Averaged clothing index (H1 vs H7, C1 vs C7).

This is probably because each occupant first attempts to adjust the thermal sensation with his or her clothes and then uses the remote control. Extremely hot sensitive or cold sensitive workers use remote controls, but many moderate workers would tolerate only the adjustment of their clothing. Thus, the setpoint temperature is determined according to the extreme occupants, which increases energy consumption. Therefore, providing a temperature setpoint with completely free control to the occupant may not always automatically optimize energy and thermal comfort.

Stop VRF in the interior zone (Case H8, C8)

If the VRF in the interior zone is stopped, and only the VRF in the perimeter zone is activated, the energy consumption is reduced for both heating and cooling. However, the impacts of heating and cooling on the dissatisfaction rate differed.

One-week temperature changes in the interior and perimeter zones for cases H8 and C8 are shown in Figure 25. The black line represents the temperature in the perimeter zone with the VRF turned on and the gray line represents the temperature in the interior zone with the VRF turned off. In the case of heating, the temperature in the interior zone is slightly low at the start of air conditioning but generally remains the same as that of the perimeter zone because of the reduced heating load. This reduced heating load can be attributed to the weak impact of the heat transfer from the exterior walls on the interior zone. Conversely, in the case of cooling, the temperature of the interior zone is ∼3°C higher than that of the perimeter zone because the generated internal heat cannot be removed via inter-zone ventilation within the perimeter zone. As a result, the thinning operation during cooling considerably reduces thermal comfort.

Figure 25. Temperature changes in the interior and perimeter zones (H8 and C8).

Discussion

The development of an emulator to assess the skill of operating VRF systems was reported in this study.

The emulator combined a building heat load calculation model, VRF model, and office worker model, which allowed us to evaluate the operation in terms of both energy and comfort. This was not feasible because many of the proposed VRF system models concurrently failed to evaluate human thermal comfort. Thus, the evaluation of VRF systems in terms of both energy performance and comfort was realized. However, the proposed emulator has certain limitations, as shown below.

First, the model is fairly simple in terms of how personal controllability, the most significant feature of VRF systems, affects comfort. It does not explain the effect on the psychology of office workers when the control of fan speed or air flow direction as well as setpoint temperature is allowed. As explained in Section 5.1, further research is required before such a comprehensive model becomes available.

Second, the airflow is not represented in detail. The model used in this study cannot represent the interference of air supply from multiple indoor units. Further, it also cannot represent the air supply reaching zones other than the zone where the indoor units are installed. These phenomena are highly likely to occur in reality and affect the energy performance. In addition, because the airflow distribution within a zone is not known in detail, draft calculations are also modeled in an approximate manner, as explained in Section 5.2. One pragmatic approach would be to perform a numerous CFD calculations in advance and maintain them as a table, which can then be referenced. Although the generality of the model will be lost, the calculation speed and detail can be increased.

Third, the strength of the dissatisfaction is not reflected in the score. The emulator counts the same one unit of dissatisfaction whether it is a small dissatisfaction that deviates slightly from the optimal thermal environment for the occupants or a large dissatisfaction that can lead to illness. Thus, at night or on holidays when the number of occupants is low, a high score may be achieved by ignoring them and not air-conditioning the office at all. This is not ethically acceptable. One solution would be to use the deviation from the optimal thermal environment for each occupant as a rating scale. However, unlike dissatisfaction, in the real world, the measurement of such a concept is difficult.

Systems that combine buildings, machines, and people are complex, and a single operational change can affect several aspects of the system. In many cases, there is a tradeoff between reducing energy consumption and improving comfort, and it is difficult to simultaneously improve both at the same time. There are also trade-offs between factors that bring about comfort, as suggested in the case studies in this study; for example, increasing fan speed increases the risk of drafts in exchange for a more stable room temperature.

It must be recognized that there is no great value in organizing these complex interactions and defining single standard operating procedures. This is because of numerous system combinations. For example, in the case of the physical properties of moist air, because physical phenomena are invariant, once a tool called a psychrometric chart is created, there is almost no need to update it. However, the energy and comfort characteristics of a building with VRF vary with the building structure, weather conditions, and thermal preferences of the occupants. In extreme cases, the addition of a single indoor unit significantly affects the characteristics of the system and changes its optimal operation.

Therefore, a single standard operation that can be commonly applied to various VRF systems is not suitable; rather, a meta-level methodology that enables the determination of the best operation for a specific system with each change is required.

The VRF system emulator proposed in this study is intended to be used as a testbed for verifying methods to search for optimal operations. The World Championship described in this paper is a specific application of emulators. As with the first championship held in 2019, we plan to collect information from the participants after the event, organize their ideas, how they changed their operations, and what results they obtained and report them in a paper.

To use the emulator as a teaching tool, tutorial materials explaining the modification of the operation and the effect achieved would also be necessary. We expect that knowledge for this purpose will be gained through the championship, and we plan to publish the results in the future.

The source code and the manual for the emulator described in this paper is available at www.wccbo.org. The second championship will be held in the summer of 2024, and readers of this journal will be invited to participate to accumulate knowledge on how to optimize VRF systems and promote research in this field.

Nomenclature

Acc	=	acceleration rate (-)
almt	=	cut-off parameter (-)
Ca	=	heat transfer coefficient between upper and lower zones (W/(m2·K))
cpa	=	specific heat of air at constant pressure (J/(kg·K))
D	=	number of dissatisfied occupants (persons)
DR	=	averaged total dissatisfied rate of occupant (-)
DRR	=	dissatisfaction reduction ratio (-)
E	=	energy consumption (MJ)
ECI	=	energy saving and comfort index (-)
ERR	=	energy reduction ratio (-)
E2R2	=	effective energy reduction ratio (-)
Ho	=	width of slot diffuser (m)
h	=	heat (W)
Kp	=	centerline velocity constant dependent on outlet type and discharge pattern (-)
P	=	percentile (%)
Rht	=	heat distribution ratio to the lower zone (-)
Roc	=	occupancy ratio (-)
r	=	radius (m)
s	=	length of jet stream (m)
sth	=	thermal sensation (-)
Turb	=	power of turbulence (-)
t	=	temperature (°C)
u	=	velocity of jet stream (m/s)
V	=	volumetric air flow (m3/s)
w	=	weigh factor (-)
x	=	horizontal distance from diffuser (m)
y	=	vertical distance from diffuser (m)
αt	=	thermal diffusivity (m2/h)
δ	=	standard deviation of velocity (m/s)
θ	=	angle between um and x-axis (radian)
γ	=	density (kg/m3)

Subscript

a	=	air
an	=	annulus
DR	=	dissatisfied rate
drft	=	draft
e	=	environment
E	=	energy consumption
jet	=	jet stream
L	=	Lower zone
lcl	=	local
m	=	middle (centerline)
r	=	reference
th	=	thermally
U	=	Upper zone
zn	=	zone

Acknowledgements

This research was conducted as part of the activities of the ‘BEMS Committee’ at ‘The Society of Heating, Air-Conditioning and Sanitary Engineers of Japan’.

Data availability statement

The data that support the findings of this study are openly available in [shizuku] at [https://github.com/et0614].

Additional information

Funding

This work was partially supported by Japan Society for the Promotion of Science [KAKENHI Grant Numbers: JP 23K04148], and results obtained from the project ‘Research and Development for Total Cost Reduction of Heat Utilization as Renewable Energy’, commissioned by New Energy and Industrial Technology Development Organization (NEDO) in Japan [Grant number: P19006].