Fabrication of High-Quality Polymer Composite Frame by a New Method of Fiber Winding Process

Polymer composite frame has been frequently used in the main structural body of vehicles in aerospace, automotive, etc., applications. Manufacturing of complex curved composite frame suffer from the lack of accurate and optimum method of winding process that lead to preparation of uniform fiber arrangement in critical location of the curved frame. This article deals with the fabrication of high-quality polymer composite frame through an optimal winding of textile fibers onto a non-bearing core frame using a fiber-processing head and an industrial robot. The number of winding layers of fibers and their winding angles are determined based on the operational load on the composite structure. Ensuring the correct winding angles and thus also the homogeneity of fibers in each winding layer can be achieved by using an industrial robot and by definition of its suitable off-line trajectory for the production cycle. Determination of an optimal off-line trajectory of the end-effector of a robot (robot-end-effector (REE)) is important especially in the case of complicated 3D shaped frames. The authors developed their own calculation procedure to determine the optimal REE trajectory in the composite manufacturing process. A mathematical model of the winding process, matrix calculus (particularly matrices of rotations and translations) and an optimization differential evolution algorithm are used during calculation of the optimal REE trajectory. Polymer composites with greater resistance to failure damage (especially against physical destruction) can be produced using the above mentioned procedure. The procedure was successfully tested in an experimental composite laboratory. Two practical examples of optimal trajectory calculation are included in the article. The described optimization algorithm of REE trajectory is completely independent of the industrial robot type and robot software tools used and can also be used in other composite manufacturing technologies.


Introduction
In the past few decades, polymer composites have increasingly replaced traditional materials such as wood, iron, and steel in advanced industrial applications, due to their superb mechanical features such as flexible design capability, high strength to weight ratio, thermal resistance, etc. [1,2]. Polymer composite frames are mainly used to reinforce the chassis, body, and doors of a car, or to strengthen the fuselage and attach windows to the fuselage, etc., in the aerospace and automotive industries, as well as in the manufacture of agricultural machinery [3][4][5][6]. These composite frames are primarily used due to their excellent mechanical and physical properties such as resistance to harsh weather conditions, as well as long term resistance to corrosion in severe environmental conditions, the wound fibers (usually carbon, aramid, or glass fibers) and the fabrication process [11][12][13]. Ensuring the correct fiber winding angles from the geometric point of view and hence the homogeneity of the windings is one of the important aspects of composite quality. Previous studies indicated that any inhomogeneity during fiber windings process in the fabrication process polymer composite frame, results in preparation of defect as a potential source of defect that induce stress concentration and early failure phenomena [14][15][16][17].
Polymer composite frames are normally fabricated in complex irregular geometry with various configurations (e.g., Figure 1), in which industrial robots play an important role in the production of the fiber winding process [15]. Once the core frame is wound by the fibers, then it is replaced in a mold where matrix material, such as resin, etc., is injected around the frame to form a solid layer with specific thickness as the polymer composite frame [15,18]. The advantages of robotic fiber over manual fiber windings (winding of fibers by production worker without the use of a robot or other textile machine) in the manufacturing process of composite frame were investigated by Shirinzadeh, et al. [19]. It is proven that determining a correct off-line REE trajectory results in making high-quality winding of fibers onto a core frame during the production process of the polymer composite frame [3,4,16,20]. The process of trajectory calculation during the winding process for simple frame geometry in the form of two-and three-dimensional (2D and 3D) cases, is described elsewhere [21]. However, this trajectory process is not optimized and may not be applicable in the case of a more complicated 3D shaped non-bearing core frame. Sofi et al. presented dry fiber winding possibilities and explanation of the most important processes of winding [22], and Polini et al. pointed out that the tension of winding during robotic filament winding technology is a very important parameter that influences directly the defects and the mechanical property of composite [23]. Azevedo et al. investigated the effects of mosaic winding pattern on carbon-fiber composite cylinder fabricated using filament winding, and found the optimum radius to thickness ratio for optimum strength and stiffness properties [24]. Many authors studied the problem of the correct winding of fibers on a non-bearing core frame. This problematic field is very topical and essential for the industrial production of composite frames. Along with specific selections of layers sequence and angle that is made through design process, the "correct winding angle process" is to ensure the arrangement of the fibers such that the angle of the fibers to the frame remains uniform, homogeneous, and consistent in circular-helix cross-sectional form throughout the complex 3D curved geometry of the frame. Determination of the correct trajectory of a robot during winding of the fibers is presented, for example, in [19,[25][26][27]. Gao et al. proposed highspeed fiber placement technology in a new methodology of motion planning in a redundant robotic system [25], in which the problem of time optimization of robot motion is solved [25,28,29]. The optimization of robot trajectory in the fiber winding process is also addressed, in which graph theory is used to obtain the optimal robot trajectory as well as special algorithm such as genetic algorithm, harmony search and also Bézier curves [28][29][30][31][32]. A similar topic, however, focused rather on the study of trajectory control of an arbitrary shape winding mandrel in 3D circular braiding is explored in [20,23,27]. Sofi et al. presented dry fiber winding possibilities and explanation of the most important processes of winding [22], and Polini et al. pointed out that the tension of winding during robotic filament winding technology is a very important parameter that influences directly the defects and the mechanical property of composite [23]. Azevedo et al. investigated the effects of mosaic winding pattern on carbon-fiber composite cylinder fabricated using filament winding, and found the optimum radius to thickness ratio for optimum strength and stiffness properties [24]. Many authors studied the problem of the correct winding of fibers on a non-bearing core frame. This problematic field is very topical and essential for the industrial production of composite frames. Along with specific selections of layers sequence and angle that is made through design process, the "correct winding angle process" is to ensure the arrangement of the fibers such that the angle of the fibers to the frame remains uniform, homogeneous, and consistent in circular-helix cross-sectional form throughout the complex 3D curved geometry of the frame. Determination of the correct trajectory of a robot during winding of the fibers is presented, for example, in [19,[25][26][27]. Gao et al. proposed high-speed fiber placement technology in a new methodology of motion planning in a redundant robotic system [25], in which the problem of time optimization of robot motion is solved [25,28,29]. The optimization of robot trajectory in the fiber winding process is also addressed, in which graph theory is used to obtain the optimal robot trajectory as well as special algorithm such as genetic algorithm, harmony search and also Bézier curves [28][29][30][31][32]. A similar topic, however, focused rather on the study of trajectory control of an arbitrary shape winding mandrel in 3D circular braiding is explored in [20,23,27].
The present study focuses on the production of high-quality polymer composite frame by calculation of the correct winding angles. In this respect, the winding technology, individual components used in the winding process and the mathematical model of the winding process are described in detail to address the use of a mathematical model, matrix calculus and optimization process. A 3D complex frame is considered for calculating the optimal 3D REE trajectory, as shown in Figure 1. The attempt is to find the correct winding angles and appropriate homogeneous distribution of the fiber windings through this optimization procedure, even in the case of a highly complicated frame geometry. It should be emphasized that this study deals with the issue of defining a suitable cost function in order to find the optimal off-line REE trajectory during individual steps of the frame's passage through the fiber-processing head to ensure the correct winding of the fibers onto the frame. A differential evolution algorithm is used to obtain a minimum cost function as well as finding the optimal REE trajectory. The possibility of frame collisions with a fiber-processing head is also tested during each passage step. The procedure described in this paper allows ensuring that the geometrically correct required winding angles are maintained during the winding process.

Manufacturing of Polymer Composite Frame
Technology of polymer composite frame production is a complex process, in which the shape and function of the composite frame (open ( Figure 1) or closed) are the determining parameters for the choice of material and technology of the composite parts [33]. The choice of fibrous material depends on the required physical and mechanical properties of final product design based on specific operational load and boundary condition. It is possible to use dry fibers of carbon, glass, basalt, aramid, a combination of these fibers so-called hybrids or a combination these fibers with thermoplastic fibers. The aim of winding of dry fibers is to create directionally oriented layers of fibers so that the fibers are wound homogeneous, regularly, and evenly in each layer. The number of winding layers and the angular orientation of the individual winding fiber layers are usually determined using mathematical models and simulation of the composite frame [1,[33][34][35][36][37]. The winding process is generally done through "manual-robotic winding" in which the robot is programmed manually (teach-in method). This method is time-consuming and does not guarantee a quality winding process for more geometrically complex frame shape.
There are a few procedures for applying the fiber reinforcement to the core frame, such as methods for overlaying the fibers on the core frame, and winding the fibers from the coils while rotating the core frame [1]. Most of the fiber winding procedures are applicable to constant fiber deposition without possible optimization of the fiber laying angles in 3D. The procedure described in this article allows to optimize geometrically fiber winding on 3D open and closed frames, in which the production process is shown in Figure 2a. In the first step of the production of a composite frame, a mold ( Figure 2b) with the geometry identical to the final product is fabricated to make a non-bearing and lightweight core frame (generally made of porous polyurethane foam, brown color frame in Figure 2c). Then, the fiber is wound around the core frame (Figure 2d), and the frame is inserted into the preheated mold and then the matrix (made of thermoplastics (polyurethane, etc.) or thermosets (epoxy, polyester, etc.) are injected into the mold under controlled pressure and temperature for the curing process as dictated through the vacuum injection technology. (a) production process of long-fiber reinforced polymer composite frame, (b) a laboratory mold for holding and curing of the polymer composite frame (example of wind turbine blade), (c) the core before fiber winding, and (d) wound core with fibers, positioned in the mold before injection of the polymer matrix.

Fiber Winding Geometry
This section describes a geometric interpretation of the execution of the fiber winding on a nonbearing core frame of circular cross-section.
The right-handed Euclidean coordinate system E3 is taken into account. Vectors and matrices are written in a homogeneous form (i.e., any point ,0 u , in more detail see [38]). The Euclidean norm u of vector u , where The winding of one fiber on a composite frame of a circular cross-section is from a geometrical point of view generally a helical motion. The winding fiber creates the helix on the composite frame.
This motion is the composition of rotation of a fiber around the internal axis o of the frame and its parallel translation in direction of axis o . In this way, fiber forms a helix on the surface of the frame.
Right-handed and left-handed helices are shown in Figure 3. Points R A and L A are the initial points of helices created by the winding process. (a) production process of long-fiber reinforced polymer composite frame, (b) a laboratory mold for holding and curing of the polymer composite frame (example of wind turbine blade), (c) the core before fiber winding, and (d) wound core with fibers, positioned in the mold before injection of the polymer matrix.

Fiber Winding Geometry
This section describes a geometric interpretation of the execution of the fiber winding on a non-bearing core frame of circular cross-section.
The right-handed Euclidean coordinate system E 3 is taken into account. Vectors and matrices are written in a homogeneous form (i.e., any point V = [x V , y V , z V , 1] T and any vector u = (x u , y u , z u , 0) T , in more detail see [38]). The Euclidean norm u of vector u, where u = x 2 u + y 2 u + z 2 u is used. The winding of one fiber on a composite frame of a circular cross-section is from a geometrical point of view generally a helical motion. The winding fiber creates the helix on the composite frame. This motion is the composition of rotation of a fiber around the internal axis o of the frame and its parallel translation in direction of axis o. In this way, fiber forms a helix on the surface of the frame. Right-handed and left-handed helices are shown in Figure 3. Points A R and A L are the initial points of helices created by the winding process.

Fiber Winding Geometry
This section describes a geometric interpretation of the execution of the fiber winding on a nonbearing core frame of circular cross-section.
The right-handed Euclidean coordinate system E3 is taken into account. Vectors and matrices are written in a homogeneous form (i.e., any point ,0 u , in more detail see [38]). The Euclidean norm u of vector u , where The winding of one fiber on a composite frame of a circular cross-section is from a geometrical point of view generally a helical motion. The winding fiber creates the helix on the composite frame.
This motion is the composition of rotation of a fiber around the internal axis o of the frame and its parallel translation in direction of axis o . In this way, fiber forms a helix on the surface of the frame.
Right-handed and left-handed helices are shown in Figure 3. Points R A and L A are the initial points of helices created by the winding process. Let us consider the right-handed helix p R with axis o identical to axis z of the right-handed Euclidean coordinate system E 3 (Figure 3 (left)) and its parameters reduced pitch v 0 (length of translation during rotation of fiber by one radian), helix radius r (radius of frame), angle δ of slope of the helix p R defined by relation tgδ = v 0 /r (the detailed description of helix parameters can be found elsewhere [39,40]). Then the parametric equation of helix p R can be expressed in the form of p R (t) = (r cos t, r sin t, v 0 t, 1), t ∈< 0, ∞) Point A R = (r, 0, 0, 1) T is then the initial point of right-handed helix p R . The equation of the right-handed helix p R can also be expressed as the rotation of point A R = (r, 0, 0, 1) T around axis z (we suppose o ≡ z) and its translation in a positive direction of axis z [39]: In the case of one turn of the right-handed helix p R as shown in Figure 3 (left), parameter t lies within interval t ∈< 0, 2π >. Please note that the positive winding angle ω of one fiber on composite frame is (Figure 4 and if the value of reduced pitch v 0 is close to ∞, then ω = 0 (in this case the fiber is laid on the frame parallel with internal axis o ≡ z of the frame).  The parametric equation of left-handed helix L p can be expressed as and similarly by matrix calculus as in the case of the right-handed helix in relation (1).
In the production of the polymer composite frame, three layers of fibers during the winding process and twelve strands of fibers are usually wound on the frame in each layer. One given winding layer with positive winding angle ω is considered. Then during the winding process each strand (of the total number of twelve strands) creates right-handed helix R p on the composite frame with positive angle δ of slope of the helix. In accordance with relation (2), A circle with its center in origin O of right-handed Euclidean coordinate system E3 and radius In the case of winding negative angle ω of one fiber on the composite frame, left-handed helix p L is considered and angle ω is (Figure 4 (right)) The parametric equation of left-handed helix p L can be expressed as p L (t) = (r cos t, −r sin t, v 0 t, 1), t ∈< 0, ∞) and similarly by matrix calculus as in the case of the right-handed helix in relation (1).
In the production of the polymer composite frame, three layers of fibers during the winding process and twelve strands of fibers are usually wound on the frame in each layer. One given winding layer with positive winding angle ω is considered. Then during the winding process each strand (of the total number of twelve strands) creates right-handed helix p R on the composite frame with positive angle δ of slope of the helix. In accordance with relation (2), ω = π 2 − δ. Let angles λ i are determined by relation λ i = (i − 1) 2π/12 for i = 1, 2, . . . . , 12.
A circle with its center in origin O of right-handed Euclidean coordinate system E 3 and radius r is considered. Then arms of oriented angles λ i create on the circle the vertices A i = (r cos λ i , r sin λ i , 0, 1) T of a regular twelve-rectangle ( Figure 5 (left)). and similarly by matrix calculus as in the case of the right-handed helix in relation (1).
In the production of the polymer composite frame, three layers of fibers during the winding process and twelve strands of fibers are usually wound on the frame in each layer. One given winding layer with positive winding angle ω is considered. Then during the winding process each strand (of the total number of twelve strands) creates right-handed helix R p on the composite frame with positive angle δ of slope of the helix. In accordance with relation (2) Helix p Ri (t) = (r cos( t + λ i ), r sin(t + λ i ), v 0 t, 1) T can be obtained by this way of expression. If the considerations are generalized and the rotation of points A i around axis z of the coordinate system is considered, (Figure 5 (right)) with angle ψ, the expression of helix p Ri can be written in the form Individual helices p Ri of the winding layer have the same winding angle ω, but the next helix is rotated by angle π 6 relative to the previous one. The generalized parametric equation of the left-handed winding helix in the form (4) can be derived analogously.
All layers of windings (with positive, negative or zero winding angle ω) are progressively realized along the entire circumference of the non-bearing core frame in the case of the composite frame. In practical cases the composite frame can be open or closed.

Mathematical Model of Winding Process
A brief description of the mathematical model of the winding process and defined designations and abbreviations are given in this section The actual process of winding the fibers on a non-bearing frame is implemented by the fiber-processing head ( Figure 6 (right)) and by an industrial robot (in this article industrial robot KUKA KR 16-2).
Polymers 2020, 12, 1037 8 of 30 Figure 6. Different views of robot KR 16-2 with non-bearing core frame and fiber-processing head with three guide lines.
The local right-handed Euclidean coordinate system E3 of the REE ( LCS ) is also taken into account. The position of LCS toward BCS defines position of REE in BCS .
The points and vectors with coordinates in BCS and LCS are labeled with the subscript BCS and LCS in the following text, respectively.

Fiber-Processing Head
The visual presentation of the fiber-processing head ( Figure 6) in the mathematical model is given in Figure 7 (left). The coordinates of the head individual components are defined in the BCS . The first outer rotating guide line with fiber spools forms the first winding layer of fibers (usually under the angle 45°, see Figure 4 (left)) and is presented by circle The second guide line creates a layer of longitudinally laid fibers (winding under angle 0°-parallel to the frame axes  The non-bearing core frame is firmly attached to the REE, and the fiber-processing head is fixed in the workspace of the robot ( Figure 6 (left)). In the described mathematical model, the basic right-handed Euclidean coordinate system E 3 of the robot (BCS) is taken into account. This system is often called the "robot coordinate system" for industrial robots (Figure 7 (right)). Individual parts of the mathematical winding model are described in BCS. The points and vectors with coordinates in BCS and LCS are labeled with the subscript BCS and LCS in the following text, respectively.

Fiber-Processing Head
The visual presentation of the fiber-processing head ( Figure 6) in the mathematical model is given in Figure 7 (left). The coordinates of the head individual components are defined in the BCS . The first outer rotating guide line with fiber spools forms the first winding layer of fibers (usually under the angle 45°, see Figure 4 (left)) and is presented by circle The second guide line creates a layer of longitudinally laid fibers (winding under angle 0°-parallel to the frame axes  The local right-handed Euclidean coordinate system E 3 of the REE (LCS) is also taken into account. The position of LCS toward BCS defines position of REE in BCS.
The points and vectors with coordinates in BCS and LCS are labeled with the subscript BCS and LCS in the following text, respectively.

Fiber-Processing Head
The visual presentation of the fiber-processing head ( Figure 6) in the mathematical model is given in Figure 7 (left). The coordinates of the head individual components are defined in the BCS. The first outer rotating guide line with fiber spools forms the first winding layer of fibers (usually under the angle 45 • , see Figure 4 (left)) and is presented by circle k1 with center S1 BCS = [x S1BCS , y S1BCS , z S1BCS , 1]. The second guide line creates a layer of longitudinally laid fibers (winding under angle 0 • -parallel to the frame axes o BCS ) which is not important in the mathematical model. The third guide line forms the last layer of fibers (usually under the angle-45 • ) and is presented by circle k2 with center S2 BCS = [x S2BCS , y S2BCS , z S2BCS , 1]. Circles k1 and k2 have the same radius r CIRCLE . Wounded fibers by circles (guide lines) k1 and k2 on the frame create a right and left-handed helix ( Figure 3). Points S1 BCS and S2 BCS lie on axis s BCS of the fiber-processing head (Figure 7 (left)). Axis s BCS is parallel to the coordinate axis y BCS of system BCS (s BCS //y BCS ) in the model.

Non-Bearing Core Frame
The non-bearing core frame with a circular cross-section is defined by central axis o LCS and radius r TUBE (see example of frame in Figure 8). We suppose r CIRCLE > r TUBE . Central axis o LCS is entered in LCS of the REE using a discrete set of N points B(i) LCS Figure 8). Vector b2 (i) LCS characterizes the needed rotation of the frame around axis o LCS when frame passes through fiber-processing head (detail is described elsewhere [21]). In the case of a closed frame passes through fiberprocessing head).

Robot Trajectory Optimization
Robot trajectory optimization is usually required for more geometrically 2D and 3D shaped nonbearing core frames. The procedure for optimal robot trajectory calculation is described in this section. An important assumption of the mathematical model is to ensure a constant ration of the rotational angular speed of guide line of the fiber-processing head and speed of passage of the frame through the fiber-processing head. Robot external axes can be used to control rotational angular speed of the outer guide lines k1 and k2.
The actual position of the LCS of the REE with regard to the BCS is determined by six parameters listed in the "tool-center-point" (TCP). The first three values of TCP = (x, y, z, a, b, c) specify the coordinates of the origin of the LCS in regard to the BCS (see Figure 7 on the right). The last three parameters a, b and c determine the angles of the rotations of the LCS around the axis z, y, and x with regard to the BCS.

Robot Trajectory Optimization
Robot trajectory optimization is usually required for more geometrically 2D and 3D shaped non-bearing core frames. The procedure for optimal robot trajectory calculation is described in this section.
The winding of three layers of fibers is carried out in the production of the composite frame. Three consecutive layers of fibers are often wound onto the frame at angles of 45 • , 0 • and −45 • (Figure 6). Middle longitudinal layer of fibers is fastened to the frame by the third fiber winding. Please note that in the considered model, axis s BCS of the fiber-processing head is parallel to the coordinate axis y BCS . The correct winding angle and homogeneity of the winding fibers (filaments) on the frame are assured if the frame axis o BCS passes through a fictitious winding plane at the same point as the head axis s BCS and at the same time the axis o BCS is orthogonal to the winding plane at that point. The fictitious plane of the first winding layer is designated ρ1 and the fictitious plane of the third winding layer is designated ρ2, as shown in Figures 9 and 10. The general scheme of the frame passage through the head is displayed in Figure 10. The attempt is to optimize the position of the frame when passing through the winding planes 1 ρ and 2 ρ so that its position would be as close as possible to the case of a straight rod. We are The attempt is to optimize the position of the frame when passing through the winding planes 1 ρ and 2 ρ so that its position would be as close as possible to the case of a straight rod. We are gradually looking for an optimal frame passage through the fiber-processing head for N points   The exact winding angles are ensured when frame axis o BCS (or its part) is identical to axis s BCS of fiber-processing head in BCS during winding process (i.e., the frame or its part is a straight rod). Then axis o BCS is orthogonal to the imaginary planes of the fibers winding ρ1 and ρ2 ( Figure 9) in points of intersections M1 BCS ∈ s BCS and M2 BCS ∈ s BCS . Static guide line k3 creates the middle layer of longitudinally laid fibers.
The general scheme of the frame passage through the head is displayed in Figure 10. In this case, the axes o BCS and s BCS occupy different positions in the BCS.
The attempt is to optimize the position of the frame when passing through the winding planes ρ1 and ρ2 so that its position would be as close as possible to the case of a straight rod. We are gradually looking for an optimal frame passage through the fiber-processing head for N points B(i) lCS ∈ o LCS (lying on the o-axis) in the time, when point B(i) BCS ∈ σ (where σ is orthogonal plane to axes s BCS )-we speak about the i-th step of passage of frame through the fiber-processing head, the center of fiber-processing head is point H BCS ∈ σ ( Figure 10).
Points P1 (i) BCS ∈ o BCS and P2 (i) BCS ∈ o BCS are selected to meet the following conditions. The distance of points R1(i) BCS ∈ k1 (this point lies inside k1 in the intersection with axis o BCS ) and P1 (i) BCS ∈ o BCS measured as o-arc length (distance of these points on axis o BCS , see Section 2.2.2) is equal to value S1 BCS M1 BCS . At the same time, distance points R2(i) BCS ∈ k2 (this point lies inside k2 in the intersection with axis o BCS ) and P2(i) BCS ∈ o BCS measured as o-arc length is equal to value S2 BCS M2 BCS . The prescribed angles of the first and third windings of the fibers will be more accurate, the smaller the distances P1(i) BCS M1 BCS and P2(i) BCS M2 BCS will be in the i-th step of the frame passage through the winding head. We find such a position of the REE for each i (1 ≤ i ≤ N) that the location of axis o BCS in BCS (and hence the frame) will make the distances P1(i) BCS M1 BCS and P2(i) BCS M2 BCS as small as possible (in the case that the frame or its part is a straight rod both distances are equal to zero).
Using the mathematical fiber winding model described in Section 3, matrix calculus and the optimization method, the optimal trajectory of the REE is calculated, which means the optimal passage of the frame through the winding head. In this way, we ensure the correct angle of winding individual layers.
It should be highlighted that the procedure for calculating the optimal trajectory of the REE in the i-th step of the frame passage through fiber-processing head, is described in this section.
Plane σ is orthogonal to axis s BCS of the fiber-processing head (and hence also to axis y BCS ) and passes through the center H BCS of the head (Figures 9 and 10). Point H(i) BCS is randomly chosen in a defined circle with center H BCS , circle lies in plane σ. At the same time, a pair of orthogonal vectors h1 (i) BCS and h2 (i) BCS is created through the following relation (5) where Rot(z, ϕ(i)) is an orthogonal rotation matrix around axis z at angle ϕ(i) and Rot (x, ω (i)) is an orthogonal rotation matrix around axis x at angle ω (i). Elements of matrices Rot (z, ϕ (i)) and Rot(x, ω(i)) are listed elsewhere [21,38]. In relation (5) vector h1 BCS is first rotated around the z-axis by angle ϕ(i) and then around the x-axis by angle ω (i) and vector h1 (i) BCS is obtained. The sizes of angles ϕ (i) and ω (i) are randomly chosen and are within the limits set. Analogously vector h2 (i) BCS is constructed in relation (5).
In this way we randomly chose point H(i) BCS ∈ σ and mutually orthogonal vectors h1 (i) BCS and h2 (i) BCS . Then there is an unambiguously defined transformation matrix T(i) that is valid It means that matrix T(i) transforms LCS of REE to BCS that is a true relation (6). Transformation matrix T(i) is defined by relation [21,41] T(i) = L(i) · Q(i) where L(i) is translation matrix and Q(i) rotation matrix. The detailed calculation procedure of finding transform matrix T(i) is described in [21,42]. By fulfilling the conditions (6) For each i (1 ≤ i ≤ N), we seek optimal point H(i) opt BCS and vectors h1(i) opt BCS , h2 (i) opt BCS , so that distances P1(i) opt BCS M1 BCS and P2(i) opt BCS M2 BCS are the smallest possible values. It means that point H(i) opt BCS and vectors h1(i) opt BCS , h2(i) opt BCS determine the best frame position for the winding process.
Coordinates x(i) opt BCS , z(i) opt BCS define sought point H(i) opt BCS (coordinate y H(i) BCS is constant because H(i) BCS ∈ σ), sought unit vectors h1(i) opt BCS and h2(i) opt BCS are defined by angles ϕ(i) opt , ω(i) opt and by relation (5). Now, we focus on the definition of cost function F in the form where and ν 1 , ν 2 are weight constants which allow the specification of the importance of the first layer and the third layer of winding fibers. It follows from the definition of cost function F, the smaller value of function F, the better winding conditions.

Note 1
Ensuring the quality of the third winding layer is often more important than the quality of the first winding layer as it also ensures the fixing of the second placement of fibers in a longitudinal direction (second static guide line provides winding at zero angle). Therefore, constants ν 1 , ν 2 can be set ν 2 > ν 1 in the definition of cost functions F in relation (8).
We find the global minimum of cost function F defined by relation (8) in i-th step of passage frame through the fiber-processing head, i.e., Coordinates x(i) opt BCS , z(i) opt BCS define sought point H(i) opt BCS (coordinate y H(i) BCS is constant), sought unit vectors h1(i) opt BCS and h2(i) opt BCS are defined by angles ϕ(i) opt , ω(i) opt and by relation (5). As result of a calculation of x(i) opt , z(i) opt , ϕ(i) opt , ω(i) opt , transformation matrix T(i) opt is obtained. In accordance with relation (7) the transformation matrix T(i) opt is defined by relation where L(i) opt is the translation matrix and Q(i) opt is the rotation matrix of LCS toward BCS. Rotation matrix Q(i) opt in relation (10) can be decomposed in the form where Rot (z, a(i) opt ) is rotation matrix around axis z at angle a(i) opt , Rot (y, b (i) opt ) is the rotation matrix around axis y at angle b(i) opt and Rot (x, c (i) opt ) is the rotation matrix around axis x at angle c(i) opt . A detailed procedure of the calculation of Euler angles a(i) opt , b(i) opt and c(i) opt is described elsewhere [21,42].
where values x(i) opt , y(i) opt and , z(i) opt are defined by the elements of matrix L(i) opt (the detail is provided elsewhere [21,42]).

Note 2
It is necessary in the i-th optimization step to accept only such TCP (i) opt whose corresponding parameters with TCP(i − 1) opt parameters differ less than the specified limit. If the condition is not satisfied, it is necessary to seek another suitable minimum of cost function F.
Sequence TCP(i) opt , 1 ≤ i ≤ N, is calculated on the external PC and subsequently loaded into the control unit of industrial robot. In this way, the optimal REE trajectory is determined by the procedure described above. Then the control unit of robot interpolates the corresponding parameters of TCP (i) opt , 1 ≤ i ≤ N, by its internal interpolation functions similar to what described in [43]. The REE moves according to the off-line optimal trajectory thus determined.

Schematic Representation of the Procedure for Calculating the Optimal REE Trajectory
The schematic representation of calculation of sequence TCP(i) opt is described in the flowchart as shown in Figure 11. Note 3 (based on the flowchart shown in Figure 11) 1. Specification of the fiber-processing head in BCS (including coordinates of centers S1 BCS and S2 BCS of outer rotating guide lines k1 and k2 of the head, vectors h1 BCS and h2 BCS , common radius r CIRCLE of circles k1 and k2).

2.
Loading of the location of composite frame in LCS (including coordinates of points B(i) LCS , vectors b1 (i) LCS , b2 (i) LCS and values C(i) for 1 ≤ i ≤ N, radius of frame r TUBE ).

3.
Determination of more B(i) LCS points on frame axis o LCS and corresponding vectors b1 (i) LCS and b2 (i) LCS .

4.
Calculation of the optimal REE trajectory to ensure the high-quality of fiber winding on the composite frame. A differential evolution algorithm (see Section 4) is used for the optimization procedure. Determining the optimal sequence TCP(i) opt (1 ≤ i ≤ N) is the result of calculation.

5.
Storing the calculated sequence of TCP(i) opt (1 ≤ i ≤ N) in the central robot unit. Determining the optimal trajectory by linking individual corresponding parameters of consecutive following TCP(i) opt (using programming instruction of robot-linear interpolations or cubic splines).
The flowchart shown in Figure 12 describes point No. 4 of Notes 3 in more detail-the procedure for calculation of sequence ).

Note 4
The possible collisions of the composite frame and fiber-processing head (especially collisions of the frame and three guide lines k1, k3 and k2 with common radius r CIRCLE , the radius of frame is r TUBE , where r TUBE < r CIRCLE , see Figure 13) are tested in each step of passage of the frame through the head.

Use of Differential Evolution Algorithm to REE Trajectory Optimization
The procedure for finding minimum (relation (9)) of cost function F defined by relation (8) in the i-th step of frame passage through the winding head, is described. Cost function F often contains many local minima. Thus, using gradient methods for finding the global minimum of function F is not proper (the high probability that is found only local minimum). Genetic algorithm is often applied to finding optimal REE trajectory [30,44]. Hence, a differential evolution algorithm is used to minimize of cost function F . This optimization algorithm gives better results than a genetic algorithm when solving a given minimization problem. Classical differential evolution algorithm is usually denoted DE/rand/1/bin (for more detail see [45,46]). The modified differential algorithm (of DE/rand/1/bin) is used that is hereafter denoted by MDEA [46]. When using MDEA, the asymptotic convergence to the global minimum of cost function F is ensured [46,47]. The use of MDEA has already proven to be successful in optimizing in other technical areas [48]. It is often difficult to find a global minimum of cost function F in the final steps of the algorithm. However, we are then able to find a satisfactory local minimum.
MDEA is used to find

Use of Differential Evolution Algorithm to REE Trajectory Optimization
The procedure for finding minimum (relation (9)) of cost function F defined by relation (8) in the i-th step of frame passage through the winding head, is described. Cost function F often contains many local minima. Thus, using gradient methods for finding the global minimum of function F is not proper (the high probability that is found only local minimum). Genetic algorithm is often applied to finding optimal REE trajectory [30,44]. Hence, a differential evolution algorithm is used to minimize of cost function F. This optimization algorithm gives better results than a genetic algorithm when solving a given minimization problem. Classical differential evolution algorithm is usually denoted DE/rand/1/bin (for more detail see [45,46]). The modified differential algorithm (of DE/rand/1/bin) is used that is hereafter denoted by MDEA [46]. When using MDEA, the asymptotic convergence to the global minimum of cost function F is ensured [46,47]. The use of MDEA has already proven to be successful in optimizing in other technical areas [48]. It is often difficult to find a global minimum of cost function F in the final steps of the algorithm. However, we are then able to find a satisfactory local minimum.
MDEA is used to find TCP(i) opt for given value i, where 1 ≤ i ≤ N. The position of the frame in BCS is given by four parameters x H(i) BCS , z H(i) BCS , ϕ(i), ω(i). These parameters define the value of cost function F specified by relation (8). This means that every four parameters x H(i) BCS , z H(i) BCS , ϕ(i), ω(i) define one possible location of the frame in the i-th step of passage through the fiber-processing head. One individual y (i) in MDEA is defined by these four parameters and is a potential solution to this problem of finding global minimum (relation (9)) of function F. Specimen SPEC ( i ) is defined and determines the type and value ranges of each parameter of the possible individual y (i) in the i-th step of passage of the frame through the fiber-processing head. Then SPEC ( i ) can be expressed in the form [45] SPEC ( i ) = real, Lo 1 (i), Hg 1 (i) , real, Lo 2 (i), Hg 2 (i) , real 3 , Lo 3 (i), Hg 3 (i) , real, Lo 4 (i), Hg 4 (i) (12) Here denomination of real specifies that all four parameters are of a real type, values Lo j (i) and Hg j (i) specify its lower and upper limits, where 1 ≤ j ≤ 4. The specimen determines the admissible parameter values defined location of the frame in BCS.
In the MDEA we successively construct generations of individuals y(i). Each generation includes NP individuals, where each individual y(i) is a potential solution to the problem (relation (9)).
One way of forming the initial generation of individuals y(i) is given by relation y(i) m , j := Lo j (i) + rand (0, 1) · Hg j (i) − Lo j (i) (13) where index i indicates the i-th step of optimization, j determines the j-th component of the m-th individual of the initial generation (1 ≤ m ≤ NP, 1 ≤ j ≤ 4). The function rand (0, 1) randomly generates a value from a closed interval < 0, 1 >.
A sequence of generations G(k) is created, where k denotes the number of the generation. Each generation comprises of individuals y(i) and we look for the individual with the smallest value F(y(i)). In general, four individuals of the current generation of MDEA participate in the creation of an individual of the next generation. The generated individuals are saved in matrix B ∈ R NP × 5 . Each row of this matrix represents one individual y(i) and its evaluation F(y(i)). During the creation of individuals, it is essential to ensure that the components of each generated individual are consistent with relation (12).

Pseudo-Code of MDEA
In this section, detail information about MDEA is described. The Algorithm 1 consists of three basic parts: necessary input values, its own computational part, and output values of the algorithm. The algorithm is presented in pseudo-code for given fixed number i (determining the i-th step of REE trajectory), where 1 ≤ i ≤ N. The index i is not used in the following algorithm for clarity of written pseudo-code.

Input:
The number of calculated generations NG, crossover probability CR, mutation factor f , generation size NP, the dimension of individuals D = 4, lower limits Lo j and upper limits Hg j , 1 ≤ j ≤ 4. Internal computation:

1.
Create an initial generation (k = 0) of NP individuals y k m , 1 ≤ m ≤ NP, (e.g., by use of relation (13)). Find index p which satisfies the condition F y k+1 p ≥ F(y k m ) for 1 ≤ m ≤ NP, y k p := y rand , where y rand satisfies (13) end while (k).

Output:
The best found individual y opt is represented by the row of matrix B that contains the corresponding value min F y k m ; y k m ∈ B .

Comments.
The repeat until condition cycle is always executed at least once, since the controlling condition is checked at the end of the cycle. Function rand (0, 1) randomly picks a number from the interval 0, 1 . The notation y k m, j means the j-th component of an individual y k m in the k-th generation. The individual y opt in pseudo-code of MDEA is the final solution and corresponds to designation y(i) opt that includes optimized parameters x H(i) opt , z H(i) opt , ϕ(i) opt , ω(i) optn . However, it should be noted that in general parameters x H(i) opt , z H(i) opt , ϕ(i) opt , ω(i) opt calculated by MDEA can only be optimized (and not optimal) parameters in relation to equation (9). This is due to the calculation of the final number of generations of individuals using MDEA. Therefore we mark calculated parameters as optimal parameters x H(i) opt , z H(i) opt , ϕ(i) opt , ω(i) opt .

Mechanical Performance of Polymer Composite Frame
The polymer composite frames are designed and fabricated to sustain various loads under specific boundary condition. In this regard, the frame would normally face different types of tensile, compressive, shear, bending and twisting loads. Such a complex loading normally results in the development of various stresses in the composite material and cause early damage phenomena. On the other hand, the composite frames with the complex geometry in open or close form, by nature, are under stress concentration phenomena, which normally are intensified at the curved location of the frame. As mentioned in this study, one of the technological challenges in composite frame production is the homogenous fiber winding process especially at the curved section of the frame. In this regard, few images of the final winding product (before matrix injection) of an L-shape frame based on the manual-and new-robot winding processes are provided, as shown in Figure 14. Results of the new winding method indicated that the angle of the fibers to the frame remains uniform, homogeneous, and consistent in circular-helix cross-sectional at the critical curved section, while the manual-robot winding resulted in random inconsistent fiber winding at different plies. Therefore, the final product of the composite frame using manual-robot winding would contain large macro-scale sections of polymer matrix pocket (with no fiber) at various layers in the curved location where is under stress concentration. In fact, due to the superb mechanical properties of the fiber materials, the existence of matrix pocket itself, would act as a location for stress concentration, even if it was made at the straight locations of the frame that extensively weaken the composite frame [6,36,49,50].
Damage phenomena in normal composite structures (with uniform, homogeneous, and consistent fiber arrangement) under loading normally initiate with micro-scale matrix cracking/crushing that grow in size across the ply thickness (parallel to the fibers), and induce early stage of multi-delamination event. Further loading could result in local fiber cracking/buckling, excessive multiple intra-and inter-laminar failures, weaken the laminate and lead to structural rupture [13,[51][52][53][54]. The fiber arrangement in normal composites would result the fibers as the load-bearing core component of the composite, to sustain the load, and ensure the composite frame to bear the designed load before fracture. The process of damage in normal composites occurs gradually, and composite structures are normally able to sustain severe loading conditions [54,55]. However, due to the inconsistent fiber arrangement in composite frame made by manual-robot winding, and stress concentration phenomena at curved sections as well as the matrix pocket parts, the material damage that initiate as matrix micro-crack would easily shift to cross-sectional fracture of the frame at the matrix pocket locations. In this condition, due to the very low properties of matrix materials in comparison with the fiber properties (10%), it is expected that each ply could only sustain about 10% of the load-share assumed for the composite ply. This excessively diminish the yielding limit and desired tolerance of the composite frame [6,53,55,56]. the final product of the composite frame using manual-robot winding would contain large macroscale sections of polymer matrix pocket (with no fiber) at various layers in the curved location where is under stress concentration. In fact, due to the superb mechanical properties of the fiber materials, the existence of matrix pocket itself, would act as a location for stress concentration, even if it was made at the straight locations of the frame that extensively weaken the composite frame [6,36,49,50]. Damage phenomena in normal composite structures (with uniform, homogeneous, and consistent fiber arrangement) under loading normally initiate with micro-scale matrix cracking/crushing that grow in size across the ply thickness (parallel to the fibers), and induce early stage of multi-delamination event. Further loading could result in local fiber cracking/buckling, excessive multiple intra-and inter-laminar failures, weaken the laminate and lead to structural rupture [13,[51][52][53][54]. The fiber arrangement in normal composites would result the fibers as the loadbearing core component of the composite, to sustain the load, and ensure the composite frame to bear the designed load before fracture. The process of damage in normal composites occurs gradually, and composite structures are normally able to sustain severe loading conditions [54,55]. However, due to the inconsistent fiber arrangement in composite frame made by manual-robot winding, and stress concentration phenomena at curved sections as well as the matrix pocket parts, the material

Practical Experimental Verification Tests of Optimization Procedure
The calculation method of the optimal REE trajectory described in this article was applied to practical problems. If the non-bearing core frame is not too complicated a 3D shape then it is possible to use a simpler non-optimized algorithm to calculate the REE trajectory described by Martinec, et al. [21]. For more complicated 2D and 3D shaped frames it is suitable to use the algorithm for calculating of optimal REE trajectory described in this article.
Entering of the input data is simple and as the output values of the described procedure we get the resulting sequence TCP(i) opt , 1 ≤ i ≤ N. In this way, optimal REE trajectory is calculated on an external PC. Sequence TCP(i) opt is subsequently stored in the central robot unit and in this way optimal REE trajectory is defined.
The focus was on calculation of optimized REE trajectory for two non-bearing core frames with a circular cross-section. In the first experimental test, a normal 2D shape was selected, while a very complicated 3D shape was selected for the second experiment. Cost function F defined by relation (8) is used in the optimization procedure of both cases.
The three layers of the fibers are simultaneously wound on the frame at angles 45 • , 0 • and −45 • and parameters ν 1 = ν 2 = 1 in definition of cost function F as provided in relation (8).

Experimental Test 1-Composite Non-Bearing Core Frame Shaped in 2D
The considered composite frame is part of the supporting structure of a baby carriage prototype ( Figure 15).

Experimental Test 1-Composite Non-Bearing Core Frame Shaped in 2D
The considered composite frame is part of the supporting structure of a baby carriage prototype ( Figure 15). The non-bearing core frame before fiber winding is shown in Figure 16 (left). The fiberprocessing head (Figure 16 (right)) is represented in the model by circle 1 k and 2 k (they correspond to outer guide lines, see Figure 7 on the left) with centers = BCS S1 [1000, 300     The calculated optimal REE trajectory was first used for graphical simulation of the frame The calculated optimal REE trajectory was first used for graphical simulation of the frame passage through the fiber-processing head using the OfficeLite simulator software and the KUKASimPro graphic simulator of the actions of the robot (see Figure 18). Optimal REE trajectory was subsequently tested in a composite experimental workplace.           The calculated optimized trajectory lies in a plane parallel to axes z and y of BCS. Therefore, parameters x, a, and b in optimized TCP OPT are constant for the entire optimal trajectory. Figure 21 shows a graph of the optimal values F(x(i) opt , z(i) opt , ϕ(i) opt , ω(i) opt ) (see relation (9)) corresponding to the optimal trajectory REE and a graph of values F(x H BCS , z H BCS , 0, 0) corresponding to the non-optimized REE trajectory described in [21].

Experimental Test 2-3D Shape Non-Bearing Core Frame
In the second practical example, the calculation of optimal REE trajectory is performed for a 3D shaped frame shown in Figures 1 and 22 (on the left and on the right). After production, this composite frame serves to reinforce the car chassis. In the second practical example, the calculation of optimal REE trajectory is performed for a 3D shaped frame shown in Figures 1 and 22 (on the left and on the right). After production, this composite frame serves to reinforce the car chassis.    shaped frame shown in Figures 1 and 22 (on the left and on the right). After production, this composite frame serves to reinforce the car chassis.

Figures 24 and 25
show graphs of parameters of optimal TCP OPT during the winding of fibers on the frame. The graphs in Figure 24 contain the first three parameters and the graphs in Figure 25 show the last three parameters of optimal TCP OPT .
Unlike the previous experimental Test 1, the frame is 3D shaped and all six parameters of optimized TCP OPT change continuously during the winding process. The sudden changes in a and b parameters in Figure 25 do not affect the rapid changes in REE orientation. It should be noted that orientation of the REE is defined by rotation matrix Q(i) opt in relation (10). The determination of Euler angles a, b and c are not unique [21,57]. Thus, in general, different parameter values can give the same matrix Q(i) opt in the product of sub-rotation matrices (relation (11)).  Figure 24 contain the first three parameters and the graphs in Figure 25 show the last three parameters of optimal OPT TCP .       Figure 26 shows a graph of the optimal values F(x(i) opt , z(i) opt , ϕ(i) opt , ω(i) opt ) (see relation (9)) corresponding to the optimal trajectory REE and a graph of values F(x H BCS , z H BCS , 0, 0) corresponding to the non-optimized REE trajectory [21].
It follows from values of cost function F defined by relation (8) in Figure 26 (and also in Test 1- Figure 21) that the optimized calculation method for definition of REE trajectory is more suitable than the simple non-optimized method defined in [21] for complex frame shapes.
H H corresponding to the non-optimized REE trajectory [21].
It follows from values of cost function F defined by relation (8) in Figure 26 (and also in Test 1- Figure 21) that the optimized calculation method for definition of REE trajectory is more suitable than the simple non-optimized method defined in [21] for complex frame shapes.

Conclusions
Fabrication of high-quality polymer composite frame is highly dependent on the homogeneous and uniform winding process of filament fibers. In the dry fiber winding technology on a non-bearing core frame, three layers of fibers are gradually created on the frame. The quality of each fiber winding layer depends on keeping the correct winding angle and thereby ensuring the homogeneity of the wound fibers. Fulfillment of these conditions depends on determining the correct REE trajectory, especially in the case of a 3D shaped frame. The simple calculation method of an off-line REE trajectory described in [21] is suitable only for frames that do not have an excessively complicated shape.
In this study, the original calculation method of optimal REE trajectory is developed, described and experimentally tested to ensure optimal possible correct winding angles of individual fiber layers onto frame. This method provides a significant advantage over a manually entered REE trajectory (manual-robot winding). Technicians often prefer the manual setting of the REE trajectory by means of a teach pendant (control body for motions of a robot [58]). However, this approach of trajectory determination is time-consuming and in practice it is not feasible to define an optimal REE trajectory for our needs.
Using the above described optimization algorithm, is appropriate in general for any shape of frame with a circular cross-section, especially in the case of a more complicated 3D frame shape. The use of the described procedure makes it possible to ensure the desired fiber angles and hence the homogeneity of the fiber windings. Such a process makes it possible to produce a composite frame with high resistance to operational stress. As a result, the mechanical performance of composite frame

Conclusions
Fabrication of high-quality polymer composite frame is highly dependent on the homogeneous and uniform winding process of filament fibers. In the dry fiber winding technology on a non-bearing core frame, three layers of fibers are gradually created on the frame. The quality of each fiber winding layer depends on keeping the correct winding angle and thereby ensuring the homogeneity of the wound fibers. Fulfillment of these conditions depends on determining the correct REE trajectory, especially in the case of a 3D shaped frame. The simple calculation method of an off-line REE trajectory described in [21] is suitable only for frames that do not have an excessively complicated shape.
In this study, the original calculation method of optimal REE trajectory is developed, described and experimentally tested to ensure optimal possible correct winding angles of individual fiber layers onto frame. This method provides a significant advantage over a manually entered REE trajectory (manual-robot winding). Technicians often prefer the manual setting of the REE trajectory by means of a teach pendant (control body for motions of a robot [58]). However, this approach of trajectory determination is time-consuming and in practice it is not feasible to define an optimal REE trajectory for our needs.
Using the above described optimization algorithm, is appropriate in general for any shape of frame with a circular cross-section, especially in the case of a more complicated 3D frame shape. The use of the described procedure makes it possible to ensure the desired fiber angles and hence the homogeneity of the fiber windings. Such a process makes it possible to produce a composite frame with high resistance to operational stress. As a result, the mechanical performance of composite frame could be assumed based on the designed load, while the load capacity of the composite frame made of manual-robot winding is diminished drastically.
The principle of the optimization algorithm can also be applied to other manufacturing production of specific composites when it is necessary to determine the 3D trajectory of the REE.
The described optimization algorithm is completely independent of the type of the industrial robot and the robot software tool used. The optimization algorithm can be used successfully by technicians of industrial composite workplaces and also by programmers of software tools for industrial robots involved in production of composites.
At present, specialized companies offer commercial software tools to industrial robot users. These tools are developed for advanced areas of industrial production (e.g., cutting, welding, pressing, or packing). However, these tools are not usually available for special production technologies such as the manufacturing of composite frames and other types of composites [59,60]. In such cases, the described optimization algorithm for calculation of an optimal REE trajectory can be used successfully.