Control-Oriented Modelling of a 3D-Printed Soft Actuator.

A new type of soft actuator was developed by using hydrogel materials and three-dimensional (3D) printing technology, attracting the attention of researchers in the soft robotics field. Due to parametric uncertainties of such actuators, which originate in both a custom design nature of 3D printing as well as time and voltage variant characteristics of polyelectrolyte actuators, a sophisticated model to estimate their behaviour is required. This paper presents a practical modeling approach for the deflection of a 3D printed soft actuator. The suggested model is composed of electrical and mechanical dynamic models while the earlier version describes the actuator as a resistive-capacitive (RC) circuit. The latter model relates the ionic charges to the bending of an actuator. The experimental results were acquired to estimate the transfer function parameters of the developed model incorporating Takagi-Sugeno (T-S) fuzzy sets. The proposed model was successful in estimating the end-point trajectory of the actuator, especially in response to a broad range of input voltage variation. With some modifications in the electromechanical aspects of the model, the proposed modelling method can be used with other 3D printed soft actuators.


Introduction
In contrast to building robots from rigid materials in conventional robotics, soft, responsive, flexible, and compliant materials have been used to make composite materials that can mimic biological systems. The study of such responsive composites is often referred to as soft robotics, which is a specific category of the field of robotics. There has been a growing interest in these materials, which have potential applications in sensors and actuators. Functional components composed of soft and active materials can be 3D-printed. Furthermore, it was observed that monolithic structures independent of pneumatics or fluidics systems, which have their own shortcomings [1], can be developed as muscle-like actuators by means of materials such as shape memory alloys, shape memory polymers, or responsive hydrogel materials [1][2][3].
Using 3D bio-printing as part of the fabrication of soft actuators from responsive materials such as hydrogel actuators with multi-material compositions that have spatial control was recently investigated [4]. In a recent study, polyelectrolyte chitosan hydrogel was used in the 3D-printing of a soft actuator [5]. Chitosan was found to be an appropriate material for drug release, cell manipulation, and similar biological applications where the actuation and responsiveness to external stimuli is polyelectrolyte actuator is developed. Lastly, the developed model is validated via the experimental test data.

Fabrication of the Polyelectrolyte Actuator
Like other 3D-printed products, first, the actuator model was drawn in Solidworks 2016 (Dassault Systemes, Waltham, MA, USA) and then the model was imported into an EnvisionTEC GmbH Bioplotter V2.2 software program (EnvisionTEC, Gladbeck, Germany). The required materials including medium molecular weight chitosan (with 75-85% deacetylation degree) and acetic acid solution were purchased from SigmaAldrich (Sydney, NSW, Australia). A mixture of 1.6 g chitosan in 0.8 mL acetic acid solution (1 v/v%) was made under vigorous stirring at 50 • C for 2 h. The resultant 3D-print ink was then prepared through sonication and centrifugation. Lastly, the ready-to-print ink was poured into a low-temperature 3D Bioplotter syringe. For solidifying each layer after extrusion, a solution of Ethanolic Sodium Hydroxide (EtOH-NaOH) with 0.25 M NaOH (Sigma Aldrich), 70 v/v% EtOH (3:7 ratio), was prepared. The porous chitosan beam with the size of 40 mm by 8 mm by 2 mm was printed layer-by-layer (Figure 1a). The 3D printing was performed with optimized parameters of the 3D Bioplotter, as explained extensively in an earlier work by the authors [16].
of the 3D-printed polyelectrolyte actuator is developed. Lastly, the developed model is validated via the experimental test data.

Fabrication of the Polyelectrolyte Actuator
Like other 3D-printed products, first, the actuator model was drawn in Solidworks 2016 (Dassault Systemes, Waltham, MA, USA) and then the model was imported into an EnvisionTEC GmbH Bioplotter V2.2 software program (EnvisionTEC, Gladbeck, Germany). The required materials including medium molecular weight chitosan (with 75-85% deacetylation degree) and acetic acid solution were purchased from SigmaAldrich (Sydney, NSW, Australia). A mixture of 1.6 g chitosan in 0.8 mL acetic acid solution (1 v/v%) was made under vigorous stirring at 50 °C for 2 h. The resultant 3D-print ink was then prepared through sonication and centrifugation. Lastly, the ready-to-print ink was poured into a low-temperature 3D Bioplotter syringe. For solidifying each layer after extrusion, a solution of Ethanolic Sodium Hydroxide (EtOH-NaOH) with 0.25 M NaOH (Sigma Aldrich), 70 v/v% EtOH (3:7 ratio), was prepared. The porous chitosan beam with the size of 40 mm by 8 mm by 2 mm was printed layer-by-layer (Figure 1a). The 3D printing was performed with optimized parameters of the 3D Bioplotter, as explained extensively in an earlier work by the authors [16].

Electro-Chemo-Mechanical Model of the 3D-Printed Polyelectrolyte Actuator
To circumvent the complexity of multiphysics modeling, this study suggests a scalable model of the polyelectrolyte actuator. Doing so, various factors are considered in the behavior of the actuators to estimate the dynamics parameters more realistically. Thus, the actuator model comprises both electrochemical and electromechanical models considering their dynamics coupling.  Figure 2 shows an electrochemical RC model that gives the relation between applied voltage and current flow across the polyelectrolyte gel. A diffusion model to represent the current flow across the 3D-printed polyelectrolyte actuator is calculated from the Fick's law of diffusion as:

Electrochemical Modelling
where c refers to the ion concentration, y indicates the path along the thickness of the actuator, D is the diffusion coefficient constant, A is the area between the interface of electrolyte solution and

Electro-Chemo-Mechanical Model of the 3D-Printed Polyelectrolyte Actuator
To circumvent the complexity of multiphysics modeling, this study suggests a scalable model of the polyelectrolyte actuator. Doing so, various factors are considered in the behavior of the actuators to estimate the dynamics parameters more realistically. Thus, the actuator model comprises both electrochemical and electromechanical models considering their dynamics coupling. Figure 2 shows an electrochemical RC model that gives the relation between applied voltage and current flow across the polyelectrolyte gel. A diffusion model to represent the current flow across the 3D-printed polyelectrolyte actuator is calculated from the Fick's law of diffusion as: where c refers to the ion concentration, y indicates the path along the thickness of the actuator, D is the diffusion coefficient constant, A is the area between the interface of electrolyte solution and polyelectrolyte actuator, and F denotes the Faraday constant. Then, the current running across the double layer capacitance can be calculated assuming the double-layer capacitance thickness as δ:

Electrochemical Modelling
Materials 2018, 11, x FOR PEER REVIEW 4 of 13 polyelectrolyte actuator, and F denotes the Faraday constant. Then, the current running across the double layer capacitance can be calculated assuming the double-layer capacitance thickness as : . (2) Next, for a 3D-printed polyelectrolyte actuator with the thickness of the h, the diffusion equation model can be written as: Lastly, an electrochemical model relating the current flow across the actuator to the applied voltage can be obtained by simultaneously solving Equations (1)-(3) as: (4)

Electromechanical Modelling
An electromechanical model that relates the osmotic pressure caused by input voltage to the bending of 3D-printed polyelectrolyte actuator is developed in this section. Assuming the strain-tocharge ratio , the relation between the charges density and the in-plane strain in polyelectrolyte actuators can be written as [17,18]: where, for a 3D-printed soft actuator with a size of W and L as width and length, respectively, can be calculated as: Generally, actuators strains are calculated knowing the tip displacement as well as the initial length and thickness of the actuators when the bending curvature is comparatively small. However, the curvature radius is not negligible when the soft actuators experience larger bending. In this study, as shown in Figure 3 and Equation (7), the strain resulted purely from the bending of the actuator, which is obtained as: where and are the swelled length and the centered length of actuators, and refer to the bending angle and the radius of curvature of the actuator, respectively, shows the distance between the outer edge of the swelled layer and the reference plane, and is the bending curvature Next, for a 3D-printed polyelectrolyte actuator with the thickness of the h, the diffusion equation model can be written as: Lastly, an electrochemical model relating the current flow across the actuator to the applied voltage can be obtained by simultaneously solving Equations (1)-(3) as:

Electromechanical Modelling
An electromechanical model that relates the osmotic pressure caused by input voltage to the bending of 3D-printed polyelectrolyte actuator is developed in this section. Assuming the strain-to-charge ratio α, the relation between the charges density ρ ch and the in-plane strain ε in polyelectrolyte actuators can be written as [17,18]: where, for a 3D-printed soft actuator with a size of W and L as width and length, respectively, ρ ch can be calculated as: Generally, actuators strains are calculated knowing the tip displacement as well as the initial length and thickness of the actuators when the bending curvature is comparatively small. However, the curvature radius is not negligible when the soft actuators experience larger bending. In this study, as shown in Figure 3 and Equation (7), the strain ε resulted purely from the bending of the actuator, which is obtained as: where l o and l c are the swelled length and the centered length of actuators, θ and R refer to the bending angle and the radius of curvature of the actuator, respectively, y b shows the distance between the outer edge of the swelled layer and the reference plane, and K is the bending curvature of the actuator.
To estimate an accurate value of the stress induced on the bending of the actuator, the superposition of the actuations effect and the bending strain are considered together. Therefore, the total induced stress, assuming the Young's modulus of the 3D printed polyelectrolyte actuator as E, can be calculated as follows.
where the opposite signs represent the expansion and contraction of the actuator sides.
Materials 2018, 11, x FOR PEER REVIEW 5 of 13 of the actuator. To estimate an accurate value of the stress induced on the bending of the actuator, the superposition of the actuations effect and the bending strain are considered together. Therefore, the total induced stress, assuming the Young's modulus of the 3D printed polyelectrolyte actuator as , can be calculated as follows.
where the opposite signs represent the expansion and contraction of the actuator sides. Making the total moments and forces on the actuators equal zero, the bending curvature can be realized. This is somewhat related to a new defined coefficient φ as follows.
Combining Equations (9) and (6) into Equation (4), the bending curvature can be expressed as: Lastly, the model relating the applied voltage to the bending of the 3D-printed polyelectrolyte actuator can be obtained by incorporating Equation (6) and Equations (5)-(10) in the following equation [19].
Equation (11) needs to be restructured to be practical for control applications due to the presence of a hyperbolic tangent term in its denominator, which is dimensionally infinite. To deal with this, the dimensionally infinite hyperbolic tangent term can be replaced by Mittag Leffler's expansion as: where = and = 4ℎ ⁄ . This simplification confines the model in a bounded range of frequencies for input voltages within the low interest frequency range for the 3D-printed Making the total moments and forces on the actuators equal zero, the bending curvature K can be realized. This is somewhat related to a new defined coefficient ϕ as follows.
Combining Equations (9) and (6) into Equation (4), the bending curvature can be expressed as: Lastly, the model relating the applied voltage to the bending of the 3D-printed polyelectrolyte actuator can be obtained by incorporating Equation (6) and Equations (5)-(10) in the following equation [19].
Equation (11) needs to be restructured to be practical for control applications due to the presence of a hyperbolic tangent term in its denominator, which is dimensionally infinite. To deal with this, the dimensionally infinite hyperbolic tangent term can be replaced by Mittag Leffler's expansion as: where X = s and Y = D/4h 2 . This simplification confines the model in a bounded range of frequencies for input voltages within the low interest frequency range for the 3D-printed polyelectrolyte soft actuator. Hence, Equation (12) represents the model for the small values of n, so a fourth order model approximates bending displacement of the 3D-printed soft actuator as: where β is a real value and z i 's (n = 1, 2, 3) are the zeros and p i 's (n = 1, 2, 3, 4) are the poles of the soft actuator transfer function.

T-S Fuzzy System Modeling Formulation
Consider a simplified dynamic system without uncertainty systems as: Fuzzy inference rules for the T-S fuzzy dynamic model of Equation 16 can be described as follows [20].
where the matrices A i ∈ R n×n and B i ∈ R n×1 represent the subsystem parameters, z = (z 1 , . . . , z n ) T , behaviour. x is the system state variable vector and u is the system input. R i represents the i th fuzzy inference rule, Knowing M j (j = 1, . . . , n) are the fuzzy sets, the inferred fuzzy set w i can be used to calculate µ i as the normalized fuzzy membership function of inferred fuzzy sets as follows.
The singleton fuzzifier, product inference, and center-average defuzzifier are used to form the global dynamic fuzzy model of the nominal system of Equation (17) as:

Actuator Characterization
The bending of the 3D-printed polyelectrolyte hydrogel actuator can be justified by the Donnan effect phenomenon. This means that the motion of counter-ions initiated by applied voltage leads to the ionic gradient within the hydrogel networks along the direction of the electric field. This results in an osmotic pressure difference within the hydrogel structure, and, consequently, causes the deflection of a 3D-printed actuator toward the counter electrode. Several factors can be considered to characterize the behaviour of the actuator. First, the effects of the electrolyte solution and its concentration on the bending behaviour of the 3D-printed polyelectrolyte actuator should be defined. Doing so, two different electrolyte solutions, NaOH and NaCl, were used to test the 3D-printed actuator endpoint deflection behaviour. The maximum endpoint deflection for the same sizes (40 mm × 8 mm × 2 mm) and patterns of 3D-printed actuators are measured with respect to the ionic strength of the electrolyte solutions. A constant voltage of 5 V was applied between the electrodes. The concentration of NaOH and NaCl solutions ranged from 0.1 to 0.2 M with an increment of 0.2 M for each experiment. The experiments repeated three times and the results are depicted in Figure 4a. Regardless of the concentration of electrolyte solution, the actuator reached a higher deflection in the NaOH solution than in the NaCl solution. Additionally, it is observed that, for both electrolyte solutions, there were an optimum value for electrolyte ionic strength to achieve a maximum bending angle. As shown in Figure 4, this optimum value for the 3D-printed actuator tested here was near 0.12 M. From 0.1 to 0.12 M, the endpoint deflection of the 3D-printed actuator increased. This can be attributed to an increase of the free ions moving from the surrounding solution toward their counter-electrodes. However, if the concentration of the solutions were greater than the critical concentration, 0.12 M, the shielding effect of the poly-ions leads to a reduction in the electrostatic repulsion of the poly-ions, and the subsequent decrease in the endpoint deflection. Hence, in the rest of the paper, all tests are performed in NaOH electrolyte solution of 0.12 M. 0.12 M, the shielding effect of the poly-ions leads to a reduction in the electrostatic repulsion of the poly-ions, and the subsequent decrease in the endpoint deflection. Hence, in the rest of the paper, all tests are performed in NaOH electrolyte solution of 0.12 M.
Furthermore, several different patterns and dimensions of the actuators with the same material are 3D-printed (Figure 1b). Then, the behaviour of the 3D-printed actuators based on different dimensions are investigated and the results are inserted in Table 1. The results reveal that the maximum endpoint deflection of the actuators increase in proportion to their length. It also shows almost no correlation with the width of the actuator. In addition, the bending deflection was inversely correlated to the thickness of the actuator [4].
To reveal hysteresis behaviour of a 3D-printed actuator, the repeatability tests were performed. The same pattern and size of 3D-printed samples of actuators were excited with a square wave electrical stimulus. The magnitudes of the square waves were selected as 5 V with the period of 220 s and the duty cycle of 50%. The duration of the experiment was set to be 5 times the period, and the performances of the actuators were compared based on change in maximum magnitudes in each cycle, as shown in Figure 4b. The results showed that all actuators reached their maximum deflections at the first set of experiments for a specific excitation voltage.  Error bars indicate the standard deviation of the maximum distance measured over the three trials. In the next stage, the actuator's endpoint deflections were tested over time based on different input voltage magnitudes, as illustrated in Figure 5. From the results, it can be deduced that the actuation performance increased in the higher voltage. However, the occurrence of electrolyses and Maximum endpoint deflection (mm) Furthermore, several different patterns and dimensions of the actuators with the same material are 3D-printed ( Figure 1b). Then, the behaviour of the 3D-printed actuators based on different dimensions are investigated and the results are inserted in Table 1. The results reveal that the maximum endpoint deflection of the actuators increase in proportion to their length. It also shows almost no correlation with the width of the actuator. In addition, the bending deflection was inversely correlated to the thickness of the actuator [4]. To reveal hysteresis behaviour of a 3D-printed actuator, the repeatability tests were performed. The same pattern and size of 3D-printed samples of actuators were excited with a square wave electrical stimulus. The magnitudes of the square waves were selected as 5 V with the period of 220 s and the duty cycle of 50%. The duration of the experiment was set to be 5 times the period, and the performances of the actuators were compared based on change in maximum magnitudes in each cycle, as shown in Figure 4b. The results showed that all actuators reached their maximum deflections at the first set of experiments for a specific excitation voltage.
In the next stage, the actuator's endpoint deflections were tested over time based on different input voltage magnitudes, as illustrated in Figure 5. From the results, it can be deduced that the actuation performance increased in the higher voltage. However, the occurrence of electrolyses and bubbling in electrolyte caused undesirable effects and limits the application of higher voltages in such actuators. The actuator responses under various frequencies were also investigated. The square voltage of 5 V was applied between two electrodes in various frequencies including 0.0025, 0.02, 0.031, 0.054, 0.11, 0.15, and 1.1 Hz. It can be seen from Figure 6a that the maximum endpoint deflections decrease with increasing frequency since the actuators have less time to respond. In addition, response time to the first peak was measured and the results are shown in Figure 6b. From the figure, it can be observed that the actuator reaches the first peak faster as the magnitude of the peaks decreases in higher frequencies.

Modelling and Experimental Results
The Golubev method [21] is used to estimate an appropriate model for control of the actuator The actuator responses under various frequencies were also investigated. The square voltage of 5 V was applied between two electrodes in various frequencies including 0.0025, 0.02, 0.031, 0.054, 0.11, 0.15, and 1.1 Hz. It can be seen from Figure 6a that the maximum endpoint deflections decrease with increasing frequency since the actuators have less time to respond. In addition, response time to the first peak was measured and the results are shown in Figure 6b. From the figure, it can be observed that the actuator reaches the first peak faster as the magnitude of the peaks decreases in higher frequencies. The actuator responses under various frequencies were also investigated. The square voltage of 5 V was applied between two electrodes in various frequencies including 0.0025, 0.02, 0.031, 0.054, 0.11, 0.15, and 1.1 Hz. It can be seen from Figure 6a that the maximum endpoint deflections decrease with increasing frequency since the actuators have less time to respond. In addition, response time to the first peak was measured and the results are shown in Figure 6b. From the figure, it can be observed that the actuator reaches the first peak faster as the magnitude of the peaks decreases in higher frequencies.

Modelling and Experimental Results
The Golubev method [21] is used to estimate an appropriate model for control of the actuator [22][23][24][25]. Using the Golubev method for the input signal, the uncertain transfer function relating the

Modelling and Experimental Results
The Golubev method [21] is used to estimate an appropriate model for control of the actuator [22][23][24][25]. Using the Golubev method for the input signal, the uncertain transfer function relating the applied voltage to the bending angle of the actuator is estimated as: The frequency response model of the actuators with different patterns and sizes, shown in Figure 1b, is identified based on Equation (18) and the experimental data are depicted in Figure 7. It is observed that changing various parameters of the actuator, like different sizes and patterns, lead to a new set of linear voltage dependent transfer functions. Therefore, for each specific actuator, a system identification of frequency response to different voltage levels is required. Then, the T-S fuzzy model should be incorporated to interrelate the linear transfer functions at different voltage levels for improving the model accuracy.
Materials 2018, 11, x FOR PEER REVIEW 9 of 13 The frequency response model of the actuators with different patterns and sizes, shown in Figure  1b, is identified based on Equation (18) and the experimental data are depicted in Figure 7. It is observed that changing various parameters of the actuator, like different sizes and patterns, lead to a new set of linear voltage dependent transfer functions. Therefore, for each specific actuator, a system identification of frequency response to different voltage levels is required. Then, the T-S fuzzy model should be incorporated to interrelate the linear transfer functions at different voltage levels for improving the model accuracy. To identify transfer function parameters for the arbitrary pattern actuator S1, lower and upper limit input signals, 2 V (shown in Figure 8a) and 8 V are applied to the system and the range of parameters identified are as follows. Then, a two-rule based T-S model is defined as: Lastly, the outputs of the T-S fuzzy model (shown in Figure 7b) can be calculated as: A comparison of experimental tests with the T-S fuzzy model and estimated specific voltage models for 2 and 8 V input signals is shown in Figures 9 and 10. The data is based on the average of three experiments to confirm reproducibility. Additionally, the efficacy of the developed model in terms of scalability of the 3D-printed soft actuators with different patterns and sizes are shown in Figure 10b,c where two actuators with arbitrary and lattice patterns and sizes S1 and S4 are compared in response to an analogous input. These figures show supremacy of actuator end-point position estimation by T-S fuzzy modelling compared to specific constant voltage models when changing the 3D-printed actuator parameters such as pattern and size. Furthermore, looking at Figure 10b   To identify transfer function parameters for the arbitrary pattern actuator S1, lower and upper limit input signals, 2 V (shown in Figure 8a) and 8 V are applied to the system and the range of parameters identified are as follows. Then, a two-rule based T-S model is defined as: Lastly, the outputs of the T-S fuzzy model (shown in Figure 7b) can be calculated as: where µ 1 (z(t)) + µ 2 (z(t)) = 1.
A comparison of experimental tests with the T-S fuzzy model and estimated specific voltage models for 2 and 8 V input signals is shown in Figures 9 and 10. The data is based on the average of three experiments to confirm reproducibility. Additionally, the efficacy of the developed model in terms of scalability of the 3D-printed soft actuators with different patterns and sizes are shown in Figure 10b,c where two actuators with arbitrary and lattice patterns and sizes S1 and S4 are compared in response to an analogous input. These figures show supremacy of actuator end-point position estimation by T-S fuzzy modelling compared to specific constant voltage models when changing the 3D-printed actuator parameters such as pattern and size. Furthermore, looking at Figure 10b in detail reveals some discrepancy of the T-S model from experimental results, especially over longer experimental tests. This is attributed to the time-varying intrinsic nature of the polyelectrolyte soft actuator that demands the feedback controller for a compensation purpose in future study.

Conclusions
A control-oriented modelling approach for 3D-printed polyelectrolyte soft actuators was presented in this study. First, a 3D-printed actuator with an arbitrary pattern was developed and characterized based on different sizes, electrolyte concentration, input magnitudes, and frequencies. Then, a linear transfer function of the 3D-printed polyelectrolyte soft actuator was developed to estimate the actuator behavior at different voltage signals. The T-S fuzzy model was further employed for a better presentation of the actuator model in a range of voltage variations by interrelating the voltage-dependent models. The experimental results showed improved performance obtained by using the T-S fuzzy model when compared to the linear transfer function at different voltages. The proposed model could be used for other 3D-printed soft actuators with custom design geometries due to its scalability. Endpoint deflection (mm) Figure 10. (a) Input voltage signal, (b) endpoint deflection of the arbitrary pattern actuator with size of S1, and (c) endpoint deflection of the lattice pattern actuator with size of S4.

Conclusions
A control-oriented modelling approach for 3D-printed polyelectrolyte soft actuators was presented in this study. First, a 3D-printed actuator with an arbitrary pattern was developed and characterized based on different sizes, electrolyte concentration, input magnitudes, and frequencies. Then, a linear transfer function of the 3D-printed polyelectrolyte soft actuator was developed to estimate the actuator behavior at different voltage signals. The T-S fuzzy model was further employed for a better presentation of the actuator model in a range of voltage variations by interrelating the voltage-dependent models. The experimental results showed improved performance obtained by using the T-S fuzzy model when compared to the linear transfer function at different voltages. The proposed model could be used for other 3D-printed soft actuators with custom design geometries due to its scalability.