#### Abstract

Modeling for robotic joints is actually complex and may lead to wrong Pareto-optimal solutions. Hence, this paper develops a new hybrid approach for multiobjective optimization design of a flexure elbow joint. The joint is designed for the upper-limb assistive device for physically disable people. The optimization problem considers three design variables and two objective functions. An efficient hybrid optimization approach of central composite design (CDD), finite element method (FEM), Kigring metamodel, and multiobjective genetic algorithm (MOGA) is developed. The CDD is used to establish the number of numerical experiments. The FEM is developed to retrieve the strain energy and the reaction torque of joint. And then, the Kigring metamodel is used as a black-box to find the pseudoobjective functions. Based on pseudoobjective functions, the MOGA is applied to find the optimal solutions. Traditionally, an evolutionary optimization algorithm can only find one Pareto front. However, the proposed approach can generate 6 Pareto-optimal solutions, as near optimal candidates, which provides a good decision-maker. Based on the user’s real-work problem, one of the best optimal solutions is chosen. The results found that the optimal strain energy is about 0.0033 mJ and the optimal torque is approximately 588.94 Nm. Analysis of variance is performed to identify the significant contribution of design variables. The sensitivity analysis is then carried out to determine the effect degree of each parameter on the responses. The predictions are in a good agreement with validations. It confirms that the proposed hybrid optimization approach has an effectiveness to solve for complex optimization problems.

#### 1. Introduction

Along with a modern society, human people have been facing a fast increase in stroke or accidence. Therefore, robotics has received a great interest of researchers from academics and industry. If a person is subjected to the stroke, the movement’s function of arm muscles is limited. To support the disabled people, robotic systems are designed and commercialized to assist the upper limb. In general, physicians use physiotherapy to facilitate rehabilitation process. At hospitals, doctors utilize the robots.

In the state of art of rehabilitation process, assist robots and rehabilitation devices have been designed and commercialized. Robotic devices for upper-limb rehabilitation were proposed for shoulder exercises [1]. Mechanical structures and control strategies for exoskeletons for upper-limb exoskeletons were reviewed [2]. Bilateral robots for upper-limb stroke rehabilitation were studied [3]. A whole arm wearable robotic exoskeleton is used for rehabilitation and to assist upper limb [4]. A gravity-balanced exoskeleton for active rehabilitation training of upper limb was developed [5]. Developments of active upper-limb exoskeleton robots are reviewed [6]. Controller was designed for the angular trajectory of human shoulder, elbow, and wrist joints [7]. A passive upper extremity orthosis for amyoplasia was investigated [8]. Bioinspired lower extremity exoskeleton robot was designed for supporting body gait [9]. Even though existing devices were designed well, they still had shortages such as a heavy weight and large cost. The reason is because the devices must take a motor to generate a moment, a gear pairs to transfer motions, a coil spring to store and release an elastic energy, and submechanical elements. Their mechanical elements are assembled based on kinematic joints, and this results in undesired clearances. Due to the clearances, the rehabilitation systems may be subject to undesired vibrations that can cause injury back to the disable forearm. On the other hand, because of a large size and clearance, the existing devices have a complex manufacture and control and expensive cost.

Nowadays, the patients need a flexible rehabilitation process at hospital, at therapy center, or at home. Some patients with disable muscle in upper-limb desire reasonable commercial devices for training at home. Motivated from that requirement, a new flexure elbow joint is designed in this work instead of utilizing coil spring. The joint can store and release the strain energy similar to the function of traditional springs. The compliant spring is designed based on concept of compliant mechanism due to a light weight, monolithic manufacture, a low cost, and high positioning accuracy [10–14]. The proposed joint can be made from a monolithic structure by using 3D printer or wire electrical discharged machining [15–18].

In order to meet the practice requirements, the proposed joint should have a large strain energy and a wide torque being suitable for various loads. A multiobjective optimization design is preferred in this regarding [19, 20]. In general, mathematical equations are formed, and then an evolutionary optimization algorithm is applied to seek the best solutions. However, if the established mathematical equations are wrong, the predicted results are unaccurate. For these reasons, this study introduces a new data-driven multiobjective optimization technique for optimizing the performances of the proposed joint to decrease the modeling errors. A hybrid integration includes of finite element method, Kigring metamodel, and multiobjective genetic algorithm. Traditionally, most of multiobjective evolutionary algorithms (MOEAs) often give a Pareto-optimal front [21]. Unlike previous MOEAs, the proposed hybrid algorithm can generate a lot of Pareto-optimal fronts, which provides many optimal candidates, and then a best solution for real-work purpose is chosen.

The purposes of this paper are to propose a multiobjective optimal design strategy for the flexure elbow joint in terms of good static characteristics. To improve the static performances, a hybrid optimization algorithm is developed. A validation is performed to evaluate predicted results and efficiency of the proposed approach.

#### 2. Description of Mechanical Design

##### 2.1. Flexure Elbow Joint

A flexure elbow joint (FEJ) was designed to allow only one purely rotational motion. This motion was similar to the motion of actual elbow. For physically disable people after stoke, a movement of muscle is still difficult. With the purpose of support, the proposed FEJ was used as an elastic joint to store and release a strain energy during elbow rehabilitation. The FEJ was integrated with the assistive devices; the daily life activities of people after an accident or a stroke would be supported easily.

Figure 1(a) illustrates a 3D model of the proposed FEJ. It included three annular rings, separately. The ring 1 was an inner elastic element with* t*_{1}-thickness and this ring was thinnest among three rings because it could store the strain energy better. The ring 2 was a middle elastic element with* t*_{2}-thickness. The ring 3 was an outer elastic element with* t*_{3}-thickness. The thickness of ring 3 was largest in order to make a good stiffness of the joint, followed by* t*_{1} and* t*_{2}.

**(a)**

**(b)**

Figure 1(b) gives the basic dimensions of the FEJ. This joint was located at three fixed holes. In order to make a pure rotation around center O, an angle was applied. Corresponding to a rotational angle, a reaction torque* T* was appeared. Because the FEJ was an elastic element, the reaction torque was therefore proportional to the stiffness and rotational angle according to the Hook’s law as.where* T* represents the reaction torque,* K* is stiffness of FEJ, and is rotational angle.

Before conducting a computational analysis, a 3D FEM model of FEJ was built and simulated through finite element method in ANSYS software. Coarse mesh was adopted for rigid links while three compliant springs were refined and achieved a fine mesh and better analysis results. The face size method was used for meshing. The size of the selected elements was 0.3 mm. The total number of elements was 263009. The number of nodes was 464603, as seen Figure 2. To assess the quality of meshing, skewness criteria were adopted [11]. Figure 3 shows that the element metric is concentrated in region less than 0.6. It means that the quality of meshing is good.

In this analysis, a nonlinear finite element model was used to analyze the strain energy and output torque. The polyethylene material was chosen for the proposed joint. The properties of material were given in Table 1. Distributions of deformation of the rings and stress were given in Figure 4. The results showed that during the operation process, a buckling appears at each spring but it still makes a clearance of three springs without self-contact deformation among them. This guarantees a safety working operation during the assistive process. Besides, the maximum stress concentration occurred at the middle springs where red color marks, as depicted in Figure 4. Rounds and fillets for all springs were created so as to decrease the stress concentration.

As depicted in Figure 2, the FEJ was only desired to make a pure rotation around the z-axis. Its motion was replied on the elastic elements. Therefore, three rings were refined again to achieve a good accuracy of analysis. The others were assigned as coarse meshes because they no need to be computed. As given in Figure 1, there are seven design parameters of the joint, including* w*,* R*_{1},* R*_{2},* R*_{3},* t*_{1},* t*_{2}, and* t*_{3}. The factors* w*,* R*_{1},* R*_{2}, and* R*_{3} were assigned as constant values because they did not influence the responses of the proposed FEJ according to the elastic beam theory. Meanwhile, the factors* t*_{1},* t*_{2}, and* t*_{3} were considered as the design variables. According to the beam theory, each elastic element has a cross-section area of a rectangle; the stiffness of that area can be determined by the following equation [22]:where* E* is Young’s modulus of material,* w*_{e} is the width,* t*_{e} is thickness, and* l*_{e} is the length of a cross-section area of a rectangle. According to [22], the shear was assumed to be neglected in (2) and an applied force at free end was exerted on the beam.

Compliance of the FEJ was desired for elastic springs. If the compliance was increased, the strain energy was raised well. However, an increase in stiffness resulted in a decrease in compliance or strain energy. As seen in (2), the stiffness was dependent on the width, thickness, and length of elastic spring. It means that an increase in* w*_{e} and* t*_{e} leads to raise of stiffness; however an increase of the length results in a decrease in the stiffness. It could be concluded that the thicknesses* t*_{1},* t*_{2}, and* t*_{3} are the most important parameters affecting the performances of the joint. These parameters were selected as the design variables in this study. The polyethylene was selected as the material for the FEJ because its light weight. The parameters of material and design parameters were given, as in Table 1.

##### 2.2. Primary Application for Upper-Limb Assistive Device

In general, a pulley, kinematic joint, and string are used as a gravity-balanced mechanism to support movement of the forearm but costly, as illustrated in Figure 5. Instead of using traditional system, this study suggests a planar compliant springs. As shown in Figure 6, the main components of a new upper-limb assistive device are upper planar spring, lower planar spring, and flexure elbow joint. Three proposed springs make a gravity-balanced system for the forearm in which the FEJ is the most important element to generate the rotational motion for elbow. They can be fabricated by wire electrical discharged machining process as a monolithic structure. The forearm can be located by a trap while the hand approaches the handle. During rehabilitation process, the patient is located on the seat or wheelchair. A motor is used to generate a moment. This rehabilitation system can be integrated with an intelligent controller to control the assistive process. With the support from the system, a patient can pick some objects such pen, bottle, fruit, etc. The gravity-balanced system is aimed at compensating the weight of the human arm.

#### 3. Formulation of Multiobjective Optimization Problem

Commercialized devices can support well the patients but the cost is expensive. The proposed FEJ was proposed to decrease the motors and actuators. Regarding the efficiency of the system, the challenges which have been facing the FEJ includes the following: (i) a large strain energy so that it can store and release the elastic deformation and (ii) a large torque required to permit a good load capacity. These two requirements were conflicted together; a new hybrid optimization algorithm was therefore suggested to trade off them.

##### 3.1. Design Variables

It was noted that the proposed FEJ was very sensitive to the thicknesses. Therefore, the various thicknesses were considered as the design variables. The vector of design variables was set as .

The lower and upper bounds for the design variables were assigned as follows:where* t*_{i} (i =1,, 2, 3) represents the thickness of FEJ.

##### 3.2. Objective Functions

The multiple quality performances of FEJ were required as follows: (i) the strain energy, , was desired as large as possible, and (ii) the torque, , is required as high as possible. In summary, the optimization problem was briefly addressed as

##### 3.3. Constraints

The equivalent stress of the FEJ during the operation must be under the yield strength of the material to guarantee the elastic limitation, which was described as follows:where is the equivalent stress, is the yield strength of proposed material, and* n* is the safety factor. A higher safety is better for the system.

The strain energy was expected higher than 0.003 mJ while the torque was to over 500 Nmm.

Equations (7)-(8) were determined based on following mainly important criteria: according to the biomechanical engineering, assistive devices are proposed to fulfill clinical requirements so as to guarantee a good interaction between the device and patient [23, 24]. Such devices must have a lightweight and compliance for a comfort wear. The proposed FEJ is compact, compliant, and lightweight enough for the upper limb. During activities of daily living of patients, the motion’s range of the FEJ is desired as large as possible so that the patients can handle various objects, such as eating and drinking during their activities [23, 24]. If the elastic deformation of the FEJ is enhanced, the range of motion is improved. According to the theory of compliant mechanism, the elastic deformation of the FEJ is directly proportional to the strain energy [25]. With the support of the device, a large strain energy is required larger than 0.003 mJ so as to reduce the working efforts of patients during their daily living. The FEJ produces a good compliance with a large strain energy but it also requires a high output torque. The aim of producing a high torque is to enlarge the capability of the device for different patients. Elderly with weak muscles require a large strain energy and a low output torque while young patients with stronger muscles need a large strain energy and a higher output torque [23, 24]. In this study, an output torque over 500 Nmm is required to fulfill the clinical requirements of different patients and decrease working efforts of the users.

#### 4. Hybrid Optimization Algorithm

According to a basic optimization problem, mathematical models should be formed prior to implement an evolutionary algorithm. The approximate models always have errors because the models are mainly dependent on capability and knowledge of researchers about engineering and mathematics. Therefore, a solution may be wrong. To overcome this limitation, a hybrid approach of RSM, FEM, Kigring metamodel, and MOGA was proposed in this study so as improve the quality performances of the FEJ. Figure 7 illustrates a systematic flowchart for the multiobjective optimization for the proposed joint.

The optimization experienced the main steps, such as design a mechanical structure, define design variables and objective functions, build 3D model, evaluate initial performances, and establish the numerical experiments using response surface. These steps can be found in detail by [26].

After collecting the numerical data, the Kigring metamodel is used to build the relationship between inputs and outputs. There are various regression models such as full 2nd-order polynomials, artificial neural network, and nonparametric regression. In this study, Kigring metamodel was a suitable choice for the estimated database. It was considered as a black-box to approximate the complex nonlinear relationship between design parameters and the qualities. Kigring metamodel was used to find a pseudoobjective functions. These objective functions were used for the MOGA algorithm. Kigring metamodel is an interpolative Bayesian metamodeling technique combining a global model with local deviations [27]:where is a known polynomial function and is a stochastic process with mean zero and nonzero covariance. The nonzero covariance is determined:where is the correlation matrix. R() is the correlation function between any two sample points xi and xj. The Gaussian correlation function was described aswhere* n* is the number of design variables and is the unknown correlation parameter to be determined. Predicted estimates at untried values of are computed:where* y* is the column vector of length n that contains the sample data of the responses, is a column vector of length that is filled with ones when is taken as a constant, and is the correlation vector between a predicted point x and the sample points described as is estimated byThe estimated variance of output model can determined asThe unknown parameters obtained using the maximum likelihood estimation can be formulated as [27]where and are the functions of .

Accuracy of a regression model could be evaluated based on the four criteria, including the coefficient of determination, the root mean square error, the maximum relative residual, and the relative root mean square error.

The coefficient of determination (*R*^{2}) with its value in the range was determined as

Root mean squared error (*RMSE*) was defined as a measure of the differences between values predicted by either a model or an estimator and the observed actual values.

Relative maximum absolute error (*RMAR*) was defined as the difference between the observed values and the predicted values and was determined as

Relative average absolute error (*RAAE*) was determined as the difference between predicted values by a model or an estimator and the observed true values.* RAAE* measures the average of the squares of the errors.where* m* is the number of observations;* a*_{i} and* u*_{i} denote the actual and predicted, respectively.

If the regression models were not well established, this step would be refined by adjusting the range of design variables. Otherwise, it would move on the next step.

After the estimated objective functions were determined, the optimization process was implemented by programing multiobjective genetic algorithm (MOGA). This algorithm was proposed for multiobjective optimization problem in this study because it can converge to the global Pareto solutions. This algorithm helped to seek a Pareto-optimal set for multiple objective optimization [26]. The MOGA was a variant of the Nondominated Sorted Genetic Algorithm-II (NSGA-II) that replied on controlled elitism concepts. It can solve multiple objectives and constraints. At last, the MOGA can find the global optimum solution. The controllable parameters of MOGA in this study were given in Table 2. These values were chosen through many simulations.

If the optimal results were satisfied, the extra validations were conducted to evaluate the robustness and efficiency of the proposed hybrid optimization approach. The optimization process was ended herein. If they were not well found, a refinement was the most important phase to seek the optimal solutions. After optimal candidates were generated and based on the initial requirements of quality characteristics, the researchers would evaluate the candidates. If there was no any candidate that was satisfied, the ranges of quality characteristics or the range of design variables must be controlled or refined again. This step was repeated until the best candidate was found.

#### 5. Results and Discussion

##### 5.1. Data Collection

A 3D FEM model was designed and then FEA simulations were implemented to collect the performances of the proposed FEJ. The number of numerical experiments was determined using (9). The data were retrieved the FEM, a given in Table 3. And then, based on Kigring regression model, each approximated function was made for the strain energy and torque. The results indicated that the stress is still much lower than the yield strength of polyethylene. This guaranteed a long working fatigue life.

##### 5.2. Kigring Metamodels

Nowadays, there are some common surrogate models such as the nonscreening, the artificial neural network, full second order polynomials, and the Kigring metamodel. In order to establish the regression models, the fifteen sets of sample points in Table 3 were utilized. The results gave that the coefficient of determination of the Kigring model is approximately unity. This parameter was better than that of full second-order polynomials and neural network, as shown in Table 4. The rest metrics, such as the root mean square error, the relative maximum absolute error, and the relative average absolute error of the Kigring model, were smaller than those of other models. Hence, the Kigring model was adopted for this study.

##### 5.3. Contribution of Design Variable

The contribution of each parameter to the responses was analyzed by using the analysis of variance according to 95% confidence intervals. Considering the strain energy, the results showed that two thicknesses* t*_{2} and* t*_{3} are totally statistical significance with the p-value less than 0.05 expect for the thickness* t*_{1} with p-value a little higher than 0.05. It showed that the parameter* t*_{3} has a largest contribution to the strain energy with 50.78% (F-value of 22.23), followed by parameter* t*_{2} with contribution of 13.16% (F-value of 13.67), and parameter* t*_{1} has a smallest contribution of 8.84% (F-value of 4.18), as given in Table 5.

A similar way to the reaction torque**, **the results showed that three thicknesses are statistical significance with the p-value less than 0.05. Besides, the result indicated that the parameter* t*_{2} has a largest contribution to the reaction torque with 56.53% (F-value of 34.03), followed by parameter* t*_{1} with contribution of 20.96% (F-value of 12.16), and parameter* t*_{3} has a lowest contribution of 15.89% (F-value of 9.59), as given in Table 6.

##### 5.4. Optimal Results

In order to solve an optimal solution and tradeoff among between the quality responses, the MOGA was then integrated with the FEM, the RSM, and Kigring metamodels. During the optimization process, the optimal results were generated automatically. And then, if the optimized solutions are not satisfied, an extra adjustment was embedded in the proposed methodology. It was mainly dependent on the designer’s experiences. To save computational time and the further adjustment, this study proposed two ways to gain the candidates. The first assignment was to adjust the range of design parameters. The second was a change in the range of objective functions. After implementing many computational simulations, the results found that the proposed hybrid methodology could be effective only if the range of two quality responses were limited well. To suppress the dependence on the expert knowledge, two common rules were proposed as follows:(i)The lager-the better objective function, the range of this function should be overcome an allowable threshold. For example, the largest strain energy of FEJ was desired higher than 0.003 mJ, and the range of this response was therefore controlled within the upper range of 0.003 mJ, as given in Table 7.(ii)A maximum reaction torque was also required being more than 500 Nmm, and this function was herein constrained within the lower range of 500 Nmm. The optimal result could reach the desired value of larger than 500 Nmm, as shown in Table 7.

During the optimization process, the history chart of strain energy and that of reaction torque were retrieved so that they could be virtualized. As seen in Figure 8, the strain energy was convergent being close to the upper range of 0.003 mJ. As a result, the optimal strain energy could be found larger than 0.003 mJ as the design requirement. Meanwhile, the reaction torque was convergent in the upper range of 500 Nmm, as shown in Figure 9. Hence, the optimal torque could be found around 570 Nmm.

The charts were so noise because the optimal solutions are constrained in a range so as to achieve the convergence. It means that the optimal result can be found in the range of desired space. The curves may be so noise and it is true that there were so many data points occurring at the optimum space and the optimal cost functions were found in the ranges. This was a new approach to gain a convergent solution. It allowed choosing a best solution for a real-work problem.

The virtual mathematical models were found by using Kigring metamodels. To achieve the optimal results for the FEJ, the constraints for two objective functions were set up. There have been six potential candidates generated for prefer selection. Table 8 shows how to choose a best candidate and this was dependent on user’s requirements. The result revealed that the candidate 6 was chosen as the best optimal design because it fully satisfies the mentioned design objectives (see (3)-(8)). Moreover, the equivalent stress was about 20.8 MPa that was still under the yield strength with a high safety of 1.2. This could guarantee a fatigue life and long working time.

##### 5.5. Analysis of Sensitivity

Along with optimization process, a sensitivity analysis is also a necessary step to determine an influence of each design parameter on each quality response. Commonly, there are various techniques which can be applied for calculating the sensitivity as the Nelson method, Modal method, Matrix perturbation method, Differential method, and RSM. In general, the direct differential method takes more time to analyze because it needs to construct physical models. The sensitivity could be calculated by following formula: where , are the response and the design variable ith, respectively.

The RSM [28] was chosen for this analysis. Considering the local sensitivity, the thickness* t*_{3} had a highest influence or significant contribution on the strain energy, followed by the thickness* t*_{2} while the thickness* t*_{1} had a lowest contribution. It was noted that a change in* t*_{3} would adjust so as achieve the strain energy as desired. Regarding the reaction torque, the thickness* t*_{2} had a largest influence, followed by the thickness* t*_{3} and the thickness* t*_{1} had smallest contribution, as given in Figure 10.

Particularly, Figure 11 illustrates the effects of the thicknesses* t*_{1} and* t*_{2} on the strain energy and torque, as given in Figures 11(a) and 11(b), respectively. The results indicated that a decrease in the thickness* t*_{1} results in a decrease in the strain energy. But an increase in the thickness* t*_{2} leads to raise of the strain energy correspondingly. Figure 12(b) shows the effects of the thicknesses on the reaction torque. It is noted that the torque is linearly proportional corresponding to the thickness. This response was changed sharply with an increase in the thickness.

**(a)**

**(b)**

**(a)**

**(b)**

The contribution diagram of thickness* t*_{3} was plotted, as shown in Figure 12. The results revealed that a decrease in the parameter* t*_{3} leads to a decrease in the strain energy, as seen in Figure 12(a). Meanwhile, the torque was sharply increased when* t*_{3} was lowered, as given in Figure 12(b).

Summary, almost the mentioned design parameters had significant contributions on the displacement and safety factor. This would help designers and researchers to make a decision and meet the requirements of a specific system.

#### 6. FEA Verifications

To evaluate and validate the optimal results of the proposed FEJ, a few of FEM tests were carried out. The optimal design variables of Candidate 1 from Table 8 were used as the optimal values. It was used to make the 3D model (*t*_{1} = 0.97 mm,* t*_{2} = 0.9 mm,* t*_{3} = 0.51 mm). Some FEA validations were conducted in ANSYS 18 software. Materials, boundary conditions, and load were similarly in previous steps. The quality of mesh was achieved, and then the output responses were retrieved.

As given in Table 9, the results indicated that the validated results were found at the strain energy of 0.0032 mJ and reaction torque of 591.08 Nmm. The results found that the predicted and verification results are in a good agreement. It means that the proposed hybrid optimization method is an efficient approach to solve multiobjective optimization problem for the FEJ. It can be applied to solve complex optimization problems. It would support well the computational design process of mechanical element for the assistive devices.

#### 7. Conclusions

This paper has presented a new efficient multiobjective optimization approach for a flexure elbow joint. The joint was designed based on connecting in series of the leaf springs. The proposed FEJ was used to support one degree of free rotation of elbow. It was integrated into the upper-limb assistive device for disable people. The FEJ was the most important element, and the static performances was therefore optimized so as meet the requirements of the device. To improve overall static performances, including the strain energy and reaction torque, a hybrid optimization approach was developed. This approach was an integration of FEM, RSM, Kigring metamodel method, and MOGA. Six optimal candidates were retrieved and then candidate 1 was chosen as the optimum solution.

ANOVA was employed to identify the significant contribution of each variable to the responses. The sensitivity analysis through the RSM was conducted to determine influence of each factor. It found that the parameter* t*_{3} has a largest contribution to the strain energy with 50.78% (F-value of 22.23), followed by parameter* t*_{2} with contribution of 13.16% (F-value of 13.67), and parameter* t*_{1} has a smallest contribution of 8.84% (F-value of 4.18. It indicated that the parameter* t*_{2} has a largest contribution to the reaction torque with 56.53% (F-value of 34.03), followed by parameter* t*_{1} with contribution of 20.96% (F-value of 12.16), and parameter* t*_{3} has a lowest contribution of 15.89% (F-value of 9.59). The results showed that the optimal results were found at the strain energy of 0.0033 mJ and the reaction torque of 588.94 Nm. The predicted results were highly consistent with both the FEA results. It was confirmed that the proposed hybrid optimization approach is effective and reliable to solve for complex optimization engineering problems.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.03-2018.11.