Temperature and Humidity PID Controller for a Bioprinter Atmospheric Enclosure System

Bioprinting is a complex process, highly dependent on bioink properties (materials and cells) and environmental conditions (mainly temperature, humidity and CO2 concentration) during the bioprinting process. To guarantee proper cellular viability and an accurate geometry, it is mandatory to control all these factors. Despite internal factors, such as printing pressures, temperatures or speeds, being well-controlled in actual bioprinters, there is a lack in the controlling of external parameters, such as room temperature or humidity. In this sense, the objective of this work is to control the temperature and humidity of a new, atmospheric enclosure system for bioprinting. The control has been carried out with a decoupled proportional integral derivative (PID) controller that was designed, simulated and experimentally tested in order to ensure the proper operation of all its components. Finally, the PID controller can stabilize the atmospheric enclosure system temperature in 311 s and the humidity in 65 s, with an average error of 1.89% and 1.30%, respectively. In this sense, the proposed atmospheric enclosure system can reach and maintain the proper temperature and humidity values during post-printing and provide a pre-incubation environment that promotes stability, integrity and cell viability of the 3D bioprinted structures.


Introduction
Bioprinting is a booming technology that could revolutionize regenerative medicine [1]. It has a specific application in medicine of traditional additive manufacturing technology that superposes layers of material to build a biological structure. Although there are several bioprinting techniques, all of them present some common challenges to be solved, such as cell death during and after bioprinting, a long bioprinting time or insufficient micro-vascularization, among others [2]. Bioinks are commonly cell-laden hydrogels due to their good biocompatibility, so they present highly hydrated 3D networks, such as extracellular matrix, that promote oxygen and nutrient interchange [3].
It is important to know that hydrogels are very sensitive to external changes (e.g., temperature or humidity) because of their high-water content (80-90% w/v). For this reason, several processes can be associated to uncontrolled bioprinting conditions. Some of them are related to the stability and integrity of 3D bioprinting structures [4] and other ones are focused on the survival of cells during and after the bioprinting process [5]. In the first block, drying hydrogels will increase the concentration of macromolecules and will promote their crowding [3]. Furthermore, the rheological properties of hydrogels are affected by temperature, among other parameters [6]. These two problems are usually

Atmospheric Enclosure System Design
The designed atmospheric enclosure system is a sub-divided parallelepipedal with different areas: a bioprinting sub-chamber, a climatic conditions generation sub-chamber and an electronic/mechanic components sub-chamber. All designing was done using Autodesk Inventor ® , and different area dimensions are shown in Table 1. Table 1. Atmospheric enclosure dimensions.

Areas of the Atmospheric Enclosure Size (mm)
Atmospheric enclosure full size 600 × 410 × 420 Bioprinting sub-chamber (z1) 380 × 242 × 340 Climatic conditions generation sub-chamber (z2) 380 × 242 × 75 Electronic and mechanical components sub-chamber (z3) 390 × 250 × 410 Figure 1 shows the atmospheric enclosure system prototype with its atmospheric enclosure made of methacrylate. Different methacrylate widths were used in each one of the areas to control heat insulation and minimize conduction heat losses. Figure 2 shows the different sub-chambers that make up the atmospheric enclosure system.

Atmospheric Enclosure System Design
The designed atmospheric enclosure system is a sub-divided parallelepipedal with different areas: a bioprinting sub-chamber, a climatic conditions generation sub-chamber and an electronic/mechanic components sub-chamber. All designing was done using Autodesk Inventor ® , and different area dimensions are shown in Table 1.

Areas of the Atmospheric Enclosure
Size (mm) Atmospheric enclosure full size 600 × 410 × 420 Bioprinting sub-chamber (z1) 380 × 242 × 340 Climatic conditions generation sub-chamber (z2) 380 × 242 × 75 Electronic and mechanical components sub-chamber (z3) 390 × 250 × 410 Figure 1 shows the atmospheric enclosure system prototype with its atmospheric enclosure made of methacrylate. Different methacrylate widths were used in each one of the areas to control heat insulation and minimize conduction heat losses. Figure 2 shows the different sub-chambers that make up the atmospheric enclosure system.

Atmospheric Enclosure System Design
The designed atmospheric enclosure system is a sub-divided parallelepipedal with different areas: a bioprinting sub-chamber, a climatic conditions generation sub-chamber and an electronic/mechanic components sub-chamber. All designing was done using Autodesk Inventor ® , and different area dimensions are shown in Table 1.  Figure 1 shows the atmospheric enclosure system prototype with its atmospheric enclosure made of methacrylate. Different methacrylate widths were used in each one of the areas to control heat insulation and minimize conduction heat losses. Figure 2 shows the different sub-chambers that make up the atmospheric enclosure system.   Bioprinter control is divided into three different processes: bioprinter mechanical control, atmospheric enclosure control and data visualization. All three processes were individually controlled by Arduino ® (Arduino s.r.l.; 20900 Monza, Italy) UNO boards, all connected to a Raspberry ® (Rapsberry Pi Foundation; Cambridge CB2 1NF, UK) Pi 3, and the whole system was controlled by a personalized Python ® version 3.9 (Python Software Foundation; Wilmington, DE, USA) script. This script provides the user full live control over the bioprinter atmospheric enclosure parameters while temperature, humidity and CO 2 data are shown.
The control of proper atmospheric conditions inside the enclosure needs generation and sensorization of each one of the inner parameters. In this sense, temperature is generated with two 200-W electrical finned resistances, humidity is produced with cold vaporization of water in a tank using a piezoelectrical transducer and, finally, CO 2 is injected from an external AquaMedic ® (AB Aqua Medic GmbH; Bissendorf, Germany) CO 2 pressurized bottle when a Blau ® (Barcelona Marine Farm S.L.; Barcelona, Spain) (3VA, 14mA, IP 85) electro valve is opened. Likewise, the type and location of every one of the sensors used is shown in Table 2. Generation electronics is switched ON/OFF by a TONGLING (Xiamen Hongfa Electroacoustic Co; 361021 Xiamen, China) JQC-3FF-S-Z relay moduli. The main area of the atmospheric enclosure system is the climatic conditions generation sub-chamber. All generation apparatuses were placed in this area to produce the appropriate atmospheric conditions. This sub-chamber is a separated area of the enclosure connected by 4 air inlets/outlets, 2 inlets with axial fans and 2 outlets, as shown in Figure 3. The working principle of this area is to force the air to enter into the generation sub-chamber where its temperature, humidity and CO 2 are modified and then exit to the enclosure atmosphere again. This process is continuously repeated, creating an air conditioning flow until the atmospheric parameters' target values are reached. Bioprinter control is divided into three different processes: bioprinter mechanical control, atmospheric enclosure control and data visualization. All three processes were individually controlled by Arduino ® (Arduino s.r.l.; 20900 Monza, Italy) UNO boards, all connected to a Raspberry ® (Rapsberry Pi Foundation; Cambridge CB2 1NF, United Kingdom) Pi 3, and the whole system was controlled by a personalized Python ® version 3.9 (Python Software Foundation; Wilmington, Delawere, United States) script. This script provides the user full live control over the bioprinter atmospheric enclosure parameters while temperature, humidity and CO2 data are shown.
The control of proper atmospheric conditions inside the enclosure needs generation and sensorization of each one of the inner parameters. In this sense, temperature is generated with two 200-W electrical finned resistances, humidity is produced with cold vaporization of water in a tank using a piezoelectrical transducer and, finally, CO2 is injected from an external AquaMedic ® (AB Aqua Medic GmbH; Bissendorf, Germany) CO2 pressurized bottle when a Blau ® (Barcelona Marine Farm S.L.; Barcelona, Spain) (3VA, 14mA, IP 85) electro valve is opened. Likewise, the type and location of every one of the sensors used is shown in Table 2. Generation electronics is switched ON/OFF by a TONGLING (Xiamen Hongfa Electroacoustic Co; 361021 Xiamen, China) JQC-3FF-S-Z relay moduli. The main area of the atmospheric enclosure system is the climatic conditions generation subchamber. All generation apparatuses were placed in this area to produce the appropriate atmospheric conditions. This sub-chamber is a separated area of the enclosure connected by 4 air inlets/outlets, 2 inlets with axial fans and 2 outlets, as shown in Figure 3. The working principle of this area is to force the air to enter into the generation sub-chamber where its temperature, humidity and CO2 are modified and then exit to the enclosure atmosphere again. This process is continuously repeated, creating an air conditioning flow until the atmospheric parameters' target values are reached.  To ensure a sterile environment at the beginning of the bioprinting process, the bioprinter sub-chamber is sterilized by action of four UV LED panels (420 nm and 440 mW) incorporated in the bioprinter sub-chamber. After this process, air is introduced into the system from outside, passed through a HEPA filter and generates positive pressure to prevent the entry of contaminating external agents.
Although the proposed atmospheric enclosure system has a CO 2 pressurized bottle and two CO 2 sensors for future purposes, they have not been used in this study.

Mathematical Modelling
The purpose of this work is to control the inner temperature and humidity of an atmospheric enclosure system. In this sense, the inner heat balance and water mass balance equations are the general non-lineal differential equations used in temperature and humidity control processes [21,22]. In this work, terms with a negligible interaction in the total balance were not considered. Therefore, our heat and water mass balance equations [20] were given by: where C v is the volumetric air heat capacity (J/kg· • C); V is the inner air volume (m 3 ); T h is the inner temperature ( • C); G a is the air mass flow that a fan extracts from the chamber to the generation sub-chamber (m 3 /s); ρ is the air density (1.25 Kg/m 3 ); C p is the air-specific heat (1006 J/kg· • C); T c is the generation sub-chamber air temperature ( • C); Q n is the heat dissipation capacity (J/s); T 0 is the external temperature ( • C); R is the enclosure thermal resistance ( • C·s/J); d n is the inner humidity (g/kg); d c is the external humidity(g/kg); D n is the enclosure humidity gain (g/s). Transfer functions for the PID controller were obtained applying the Laplace transform to Equations (1) and (2). Defining X as 1 G a ρC p + 1 R these functions are: where G 1 (s) and G 2 (s) are temperature and humidity terms, respectively.

Experimental Evaluation
A comparative between theoretical and experimental behaviors (theoretical and experimental transfer functions) of the atmospheric enclosure was proposed. For this purpose, the different experimental behaviors of each parameter (inner temperature and humidity) were tested and determined. To obtain data for the experimental transfer function, single input, single output (SISO) and multiple inputs, multiple outputs (MIMO) tests were performed as follows: 1.
SISO test to study inner temperature when heater (electric resistance) was switched on.

2.
SISO test to study inner humidity when humidifier (cold mist humidifier) was switched on. 3.
MIMO test to study coupling of inner temperature and humidity [24,29].
Initial inputs and target values for all tests are shown in Table 3. Atmospheric values have been chosen based on criteria consulted in the literature [30]. Experimental transfer functions were obtained from test data using the Matlab ® R2018b System Identification Toolbox™. This tool is commonly used to obtain mathematical models of dynamic systems from experimental input and output data. As such, Matlab ® R2018b Simulink™ was used to simulate and compare both theoretical and experimental transfer functions.

PID Controller
Variations in temperature and humidity can change bioprinting materials' properties and affect the cellular viability, stability and integrity of 3D bioprinted structures. To minimize variations, a PID controller was added to the atmospheric enclosure system. PID controller algorithms are widely used in feedback control systems and are mathematically expressed as [31]: A good PID control system can provide the atmospheric enclosure with the necessary stability in the shortest time. In this sense, the Ziegler-Nichols closed-loop method was used in Matlab ® R2018b (MathWorks; Natick, Massachusetts, United States) to tune all PID parameters using process input/output signals. All these calculations were done through a simulation block diagram implemented in Matlab ® R2018b Simulink™, which follows the standard structure shown in Figure 4.  Experimental transfer functions were obtained from test data using the Matlab ® R2018b System Identification Toolbox™. This tool is commonly used to obtain mathematical models of dynamic systems from experimental input and output data. As such, Matlab ® R2018b Simulink™ was used to simulate and compare both theoretical and experimental transfer functions.

PID Controller
Variations in temperature and humidity can change bioprinting materials' properties and affect the cellular viability, stability and integrity of 3D bioprinted structures. To minimize variations, a PID controller was added to the atmospheric enclosure system. PID controller algorithms are widely used in feedback control systems and are mathematically expressed as [31]: A good PID control system can provide the atmospheric enclosure with the necessary stability in the shortest time. In this sense, the Ziegler-Nichols closed-loop method was used in Matlab ® R2018b (MathWorks; Natick, Massachusetts, United States) to tune all PID parameters using process input/output signals. All these calculations were done through a simulation block diagram implemented in Matlab ® R2018b Simulink™, which follows the standard structure shown in Figure  4.

Procedure
The workflow was performed as Yinping et al. [28] in their work. First, mathematical model behavior equations were solved, obtaining all constants needed to calculate the theoretical system transfer functions. Second, temperature and humidity transfer functions were calculated and evaluated with experimental transfer functions obtained from MIMO and SISO tests. Several behavior approaches of the system were obtained using the least-square method in the Matlab ® R2018b System Identification Toolbox™. All possible types in this toolbox were used: one pole (P1), two poles (P2), three poles (P3), one pole and one zero (P1Z), two poles and one zero (P2Z), three poles and one zero (P3Z), two poles sub-damped (P2U), three poles sub-damped (P3U), two poles and one zero sub-damped (P2ZU), three poles and one under-damped zero (P3ZU), one pole and one integrator (P1I), two poles and one integrator (P2I), three poles and one integrator (P3I), one pole and one delay (P1D) and two poles and one delay (P2D) transfer functions.
Next, both transfer functions' (theoretical and experimental) block diagrams were designed and evaluated. Finally, theoretical and experimental transfer functions were analyzed, obtaining the first Kp (proportional controller), Ki (integral controller) and Kd (derivative controller) values using the

Procedure
The workflow was performed as Yinping et al. [28] in their work. First, mathematical model behavior equations were solved, obtaining all constants needed to calculate the theoretical system transfer functions. Second, temperature and humidity transfer functions were calculated and evaluated with experimental transfer functions obtained from MIMO and SISO tests. Several behavior approaches of the system were obtained using the least-square method in the Matlab ® R2018b System Identification Toolbox™. All possible types in this toolbox were used: one pole (P1), two poles (P2), three poles (P3), one pole and one zero (P1Z), two poles and one zero (P2Z), three poles and one zero (P3Z), two poles sub-damped (P2U), three poles sub-damped (P3U), two poles and one zero sub-damped (P2ZU), three poles and one under-damped zero (P3ZU), one pole and one integrator (P1I), two poles and one integrator (P2I), three poles and one integrator (P3I), one pole and one delay (P1D) and two poles and one delay (P2D) transfer functions.
Next, both transfer functions' (theoretical and experimental) block diagrams were designed and evaluated. Finally, theoretical and experimental transfer functions were analyzed, obtaining the first Kp (proportional controller), Ki (integral controller) and Kd (derivative controller) values using the Ziegler-Nichols method. Then, this first PID tuning was tested and manually re-tuned to limit overshoots produced by the Ziegler-Nichols method.
To validate the atmospheric enclosure system, several tests were performed using the former calculated Kp, Ki and Kd values. Finally, to check the quality of calculated PID, an experimental error (6) was calculated.

Results
The following temperature (7) and humidity (8) theoretical transfer functions were obtained once inherent constants of the system were replaced in Equations (3) and (4).

SISO Test to Study Inner Temperature When Heater (Electric Resistance) Was Switched on
The experimental transfer function for temperature was obtained from a SISO test with laboratory environmental conditions set to 22 • C and 50% relative humidity. In the tests, actuators (electric heat resistance) were powered on or off when the inner temperature registered by sensors was below or above the target temperature (37 ± 1 • C). Data acquisition from sensors takes 700 s (test time) to assure a complete study with enough data.
The experimental transfer function, Equation (9), is calculated from the selected fitting plot for electric heater switching, the one pole and one delay (P1D) plot (Figure 5a). The accuracy rate calculated with Matlab ® R2018b for this transfer function was 99.49%.  Figure 5 shows a comparison between the experimental data (test) and the Matlab ® R2018b model from these data that extract the transfer function with a specific accuracy (model). These graphs represent the first approach of the solution that will be used as input of Section 3.3 and, after that, as initial values of the MIMO test to build a decoupled feedback of the atmospheric enclosure system.

Comparative between Theoretical and Experimental Transfer Functions
An open-loop block diagram (Figure 6) to compare theoretical and experimental SISO transfer functions before adding the PID controller was designed in Matlab ® R2018b Simulink™.

SISO Test to Study Inner Humidity When Humidifier (Cold Mist Humidifier) Was Switched on
Using previously described methodology for SISO tests, the plot that best fit the humidity behavior (Figure 5b) was one pole and one delay (P1D). Equation (10)  A slightly variation between the experimental and theoretical transfer functions (Equations (7) and (9) vs. Equations (8) and (10)) can be observed; an exponential function is multiplied to the proportional constant. Albright et al. [25] explain in their work that temperature and humidity sensors usually need this exponential function to control what they called dead time (Td), meaning a significative delay in the sensor measurement. Figure 5 shows a comparison between the experimental data (test) and the Matlab ® R2018b model from these data that extract the transfer function with a specific accuracy (model). These graphs represent the first approach of the solution that will be used as input of Section 3.3 and, after that, as initial values of the MIMO test to build a decoupled feedback of the atmospheric enclosure system.

Comparative between Theoretical and Experimental Transfer Functions
An open-loop block diagram ( Figure 6) to compare theoretical and experimental SISO transfer functions before adding the PID controller was designed in Matlab ® R2018b Simulink™. Additionally, two delay blocks have been added to this diagram to simulate the delay of sensors while reading in experimental functions. Table 3 shows all values used in inlet steps. As can be seen in Figure 7, without the addition of a controller, both temperature and humidity surpass the target values and stabilize at a higher value, taking 2500 and 250 s, respectively.  Figure 5 shows a comparison between the experimental data (test) and the Matlab ® R2018b model from these data that extract the transfer function with a specific accuracy (model). These graphs represent the first approach of the solution that will be used as input of Section 3.3 and, after that, as initial values of the MIMO test to build a decoupled feedback of the atmospheric enclosure system.

Comparative between Theoretical and Experimental Transfer Functions
An open-loop block diagram ( Figure 6) to compare theoretical and experimental SISO transfer functions before adding the PID controller was designed in Matlab ® R2018b Simulink™. Additionally, two delay blocks have been added to this diagram to simulate the delay of sensors while reading in experimental functions. Table 3 shows all values used in inlet steps. As can be seen in Figure 7, without the addition of a controller, both temperature and humidity surpass the target values and stabilize at a higher value, taking 2500 and 250 s, respectively.

MIMO Test to Study Coupling of Inner Temperature and Humidity
The dependency between temperature and humidity means that control could be a hard task when they are actuating at the same time. Wang et al. [24] exposed a higher influence of temperature on relative humidity and vice versa. They also observed a decrease by about 2-3% of relative humidity when temperature increases by 1 • C. Figure 8 shows this influence of temperature on humidity in our atmospheric enclosure system before the PID controller was set. It can be seen, without the PID controller, that both temperature and humidity fluctuate in a wide range. Micromachines 2020, 11, x FOR PEER REVIEW 9 of 16 . Figure 7. Theoretical and experimental transfer functions simulation without a PID controller.

MIMO Test to Study Coupling of Inner Temperature and Humidity
The dependency between temperature and humidity means that control could be a hard task when they are actuating at the same time. Wang et al. [24] exposed a higher influence of temperature on relative humidity and vice versa. They also observed a decrease by about 2-3% of relative humidity when temperature increases by 1 °C. Figure 8 shows this influence of temperature on humidity in our atmospheric enclosure system before the PID controller was set. It can be seen, without the PID controller, that both temperature and humidity fluctuate in a wide range. Analyzing Figure 8 plots, the system responses to some characteristic parameters of the atmospheric enclosure system are exposed in Table 4.

MIMO Test to Study Coupling of Inner Temperature and Humidity
The dependency between temperature and humidity means that control could be a hard task when they are actuating at the same time. Wang et al. [24] exposed a higher influence of temperature on relative humidity and vice versa. They also observed a decrease by about 2-3% of relative humidity when temperature increases by 1 °C. Figure 8 shows this influence of temperature on humidity in our atmospheric enclosure system before the PID controller was set. It can be seen, without the PID controller, that both temperature and humidity fluctuate in a wide range. Analyzing Figure 8 plots, the system responses to some characteristic parameters of the atmospheric enclosure system are exposed in Table 4. Analyzing Figure 8 plots, the system responses to some characteristic parameters of the atmospheric enclosure system are exposed in Table 4. Using the before-mentioned methodology, a MIMO test was performed on the atmospheric enclosure system to obtain the inherent values of the system. Laboratory environmental conditions were set to 22 • C and 50% relative humidity and the target values were the same as those used in the previous MIMO tests. Analyzing all obtained plots, P1D with delay models were the ones that best fit the temperature and humidity behaviors (Figure 9). The associated transfer functions are described in Equations (11) and (12) Using the before-mentioned methodology, a MIMO test was performed on the atmospheric enclosure system to obtain the inherent values of the system. Laboratory environmental conditions were set to 22 °C and 50% relative humidity and the target values were the same as those used in the previous MIMO tests. Analyzing all obtained plots, P1D with delay models were the ones that best fit the temperature and humidity behaviors (Figure 9). The associated transfer functions are described in Equations (11) and (12), with a 95.01% and 99.42% accuracy rate, respectively.  Figure 9 shows a comparison between the experimental data (test) and the Matlab ® R2018b model from these data that extract the transfer function with a specific accuracy (model). These graphs represent a former approach to the final solutions that will integrate other elements, such as the response time of the sensor, which will be exposed in Section 3.4.

Obtaining PID Values and Modeling System in Matlab/Simulink™
Previous results of our theoretical transfer functions cannot be considered good enough to describe our system. Interrelations between variables, non-linearity, high hysteresis or real-time variations are some reasons why temperature and humidity are complex systems to control. In this sense, performing temperature and humidity tests without any kind of control system variables will  Figure 9 shows a comparison between the experimental data (test) and the Matlab ® R2018b model from these data that extract the transfer function with a specific accuracy (model). These graphs represent a former approach to the final solutions that will integrate other elements, such as the response time of the sensor, which will be exposed in Section 3.4.

Obtaining PID Values and Modeling System in Matlab/Simulink™
Previous results of our theoretical transfer functions cannot be considered good enough to describe our system. Interrelations between variables, non-linearity, high hysteresis or real-time variations are some reasons why temperature and humidity are complex systems to control. In this sense, performing temperature and humidity tests without any kind of control system variables will provoke high fluctuation and poor stabilization of the system. Therefore, to design a good PID controller, natural fluctuations of variables were considered in the new system modelling.
Using the Matlab ® R2018b Pidtool toolbox, PID controller values of transfer functions (11) and (12) were obtained. This tool provides two different PID behavior adjustments: response velocity and long-time stabilization. PID values (Table 5) were obtained by setting a fast and robust response to the controlled system. The block diagram shown in Figure 10 uses the fitted PID values in Table 5 and compares the modelled and the real system behaviors in Matlab ® /Simulink™. The block diagram shown in Figure 10 uses the fitted PID values in Table 5 and compares the modelled and the real system behaviors in Matlab ® /Simulink™. Results from the experimental test using a decoupled PID controller of temperature and humidity can be seen in Figure 11.  Results from the experimental test using a decoupled PID controller of temperature and humidity can be seen in Figure 11. The block diagram shown in Figure 10 uses the fitted PID values in Table 5 and compares the modelled and the real system behaviors in Matlab ® /Simulink™. Results from the experimental test using a decoupled PID controller of temperature and humidity can be seen in Figure 11.  The temperature behaviors for the theoretical and the real system are very close. Figure 12 shows a maximum error for temperature of 6.01% at 185 s and the system has an average error of 1.89% once stabilized. Regarding humidity, there is a maximum error of 4.59% at 68 s, then the humidity value oscillates around the target value with an average error of 1.30% (Figure 13). The stabilization times for the theoretical models were 161 s for temperature and 281 s for humidity. Despite temperature stabilizing faster than humidity for the theoretical models, the stabilization times achieved with the proposed PID values (test values) were 311 s for temperature and 65 s for humidity. Maybe this is due to a higher error (6.01%) of the temperature model as opposed to a 1.89% error of the humidity model. Matlab ® R2018b calculates the theoretical transfer function for the temperature with a 95.01% accuracy that can justify these differences between the theoretical and experimental values.
for the theoretical models were 161 s for temperature and 281 s for humidity. Despite temperature stabilizing faster than humidity for the theoretical models, the stabilization times achieved with the proposed PID values (test values) were 311 s for temperature and 65 s for humidity. Maybe this is due to a higher error (6.01%) of the temperature model as opposed to a 1.89% error of the humidity model. Matlab ® R2018b calculates the theoretical transfer function for the temperature with a 95.01% accuracy that can justify these differences between the theoretical and experimental values.   for the theoretical models were 161 s for temperature and 281 s for humidity. Despite temperature stabilizing faster than humidity for the theoretical models, the stabilization times achieved with the proposed PID values (test values) were 311 s for temperature and 65 s for humidity. Maybe this is due to a higher error (6.01%) of the temperature model as opposed to a 1.89% error of the humidity model. Matlab ® R2018b calculates the theoretical transfer function for the temperature with a 95.01% accuracy that can justify these differences between the theoretical and experimental values.

Discussion
The objective of this work is to determine the optimum PID controller values of an atmospheric enclosure system for bioprinting. Specifically, we propose to control temperature and humidity to improve the integrity, stability and cell viability of 3D bioprinted structures. In order to design and develop our atmospheric enclosure, three components of the system were analyzed and compared with previous studies: the mathematical model, the transfer functions and the block diagram, but it was also analyzed whether to use coupled or uncoupled variables.
The similarity between the proposed atmospheric enclosure system and greenhouses or animal buildings suggests that the mathematical model used by Daskalov et al. [20] and Albright et al. [25] (heat and humidity balances) can be used in our system after some important adaptations. In this sense, the presence of animals in the study of Daskalov et al. [20] requires some factors in the mathematical model, such as sensitivity to heat or water steam produced by animals and the heat lost due to mechanical ventilation. All of these factors must be eliminated for our bioprinting environment. In addition, Albright et al. included terms of solar radiation heat, ventilation volumetric flow and plant evapotranspiration, which have no purpose in our mathematical model. Lastly, our transfer function calculations from differential equations were based on Yinping et al.'s [28] work, using a similar structure for transfer functions relating to the behavior of temperature and humidity.
According to the bibliography, there is no agreement about the best way to control temperature and humidity. So, some authors propose coupled controllers [28,32,33] while others propose decoupled controllers [24,25,34]. In this sense, knowing that an increase of 1 • C can decrease humidity by about 2-4 % [24], we performed SISO and MIMO tests to study each parameter independently and simultaneously, respectively. Considering our design and the MIMO test results (Figure 8), we finally decided to design a decoupled temperature and humidity controller that produces better results with simple PID calculations as well as increasing stability and robustness of the goal system [24].
Once the type of controller was chosen, the next step was the composition of block diagram, whose main element will be the system's plant transfer function for temperature and humidity. Previous studies proposed different types of transfer function to control temperature and humidity: some of them without delay [14,32] and others considering delay [24,34,35], which is a parameter totally dependent on the technical specifications (response time) of the electronic components [21]. Those authors that proposed no delay transfer functions try to simplify the system interactions, but their models are supposed to be less real [17,25,34]. On the other hand, the delay is commonly considered when complex transfer functions with high oscillatory responses and difficulties in stabilizing disturbances are present [25]. Hence, our system provides a better response using delay, specifically P1D type, similar to other previous studies [24,34,35].
Regarding the block diagram, three different blocks have been used: control, plant and delay. Although this diagram is quite common [25], in the same way as transfer function, delay inclusion or exclusion will modify the performance of the PID controller. So, some studies excluded the delay for temperature [17,22,28] while others considered it in their schemes for temperature and humidity [32,36]. In our case, an adapted version of Yiping et al. [28] that considers delays has been used, but other complex approaches with predictive algorithms should be considered for future improvements [36].
Finally, the global performance of our atmospheric enclosure system is pretty good. In this sense, low stabilization time and errors were obtained in the experimental tests. Other authors presented a temperature stabilization time of 2760 s for a 4 • C increment [14], 600 s for a 10 • C increment [34] or 246 s for a 5 • C increment [33]. Nevertheless, our enclosure system can stabilize a 12 • C temperature increment in 311 s. Regarding humidity, previous stabilization times were 420 s for a 10% increment [33] or 300 s for a 5% increment [34]. However, our system can stabilize a 40% humidity increment in only 65 s. It is important to note that different settings, such as heat devices, humidifiers or enclosure volume, can unfairly bias this comparison. For this reason, a specific comparative study should be performed when other atmospheric enclosure systems will be available for bioprinting. Meanwhile, we have exposed some differences among our work and the previously exposed studies for better understanding of our previous comparisons. In some cases, no technical specifications of sensors/devices are available; in this case, we cannot go into detail about the reasons of differences [34]. Other studies with a larger volume of controllers and heaters/sensors are differentially dimensioned [33]. Maybe the most similar enclosure was an incubator with comparable electronics, but unfortunately, the authors did not control humidity in the analysis [15].

Conclusions
In this work, a PID controller for temperature and humidity control of an atmospheric enclosure system for bioprinting was designed, developed and tested. Theoretical and experimental transfer functions for temperature and humidity were calculated and verified. The results show that the proposed atmospheric enclosure system is capable of stabilizing temperature and humidity in 311 s and maintaining target values with an average error of 1.89% and 1.29% for temperature and humidity, respectively. Hence, the proposed atmospheric enclosure system for bioprinting could improve post-printing environmental conditions to increase the integrity, stability and cell viability of 3D bioprinted structures.
As commented, to guarantee a proper atmospheric enclosure system for bioprinting, a CO 2 control is needed in addition to temperature and humidity. In this sense, the next step of this study will be analyzing the system's behavior after the inclusion of carbon dioxide as well as the interrelation with temperature and humidity. Additionally, bioprinters usually have a heat/cool cartridge in the extruder to control the bioink temperature. So, the addition of this heat/cool source inside the atmospheric enclosure to simulate bioprinter performance while extruding will be another future development.