Design Optimization of a Novel Networked Electromagnetic Soft Actuators System Based on Branch and Bound Algorithm

This paper presents a design optimization framework based on Branch and Bound Algorithm for novel networked electromagnetic soft actuators. The soft actuators work based on the operating principle of solenoids but are made of intrinsically soft materials. We confirmed that by scaling down the size of the soft actuators, their force to volume ratio increases. This suggests that by miniaturizing the actuators size and attaching them as a network based on the arrangement of actin and myosin filaments in skeletal muscles, the output force can be enhanced. In order to achieve the maximum output force, design parameters of a single soft actuator as well as those of a network are considered as design variables. The maximum available volume (thickness, width and length) and deflection range of the network are considered design constraints. The cost function, i.e. the output force is a non-linear mixed-integer function. A Branch and Bound optimization algorithm based on interval analysis is then proposed to solve the optimization problem. Numerical simulations are presented for a representative example of an active soft brace for the human elbow joint. The results suggest that an electromagnetic soft actuator network can provide sufficient torque to be used as a drive train for an active elbow brace for both flexion and extension over a range of around 93 degrees.


I. INTRODUCTION
In the context of physical Human-Robot Interaction pHRI, the technology requirement to fabricate the robotic platform is fundamentally different than those developed for industrial robots [1]. The fact that pHRI robots, such as industrial exoskeletons, prosthetic devices and surgical robots physically interact with the soft tissues in the human body, exerting forces outside the tissues structural limitations, introduces medical hazards. Current industrial biomedical robots are designed for fast and accurate position control applications in a perfectly controlled environment without taking into account the force interaction uncertainties that exist in pHRI The associate editor coordinating the review of this manuscript and approving it for publication was Engang Tian .
Actuators, as the source of producing force into robotic platforms, plays critical role in realizing the softness of the robot. As a result, soft actuators are gaining interests among the researchers and consequently many different types of soft actuators have been developed [5]. Soft pneumatic actuators [6]- [8] have been the first generation of soft actuators where the air pressure inside a soft and stretchable body, usually made of Polydimethylsiloxane (PDMS) [9], can be adjusted and inflate-deflate the robot's body to create a desired type of motion and amount of force. This type of soft robots are fast in response, can generate considerably high amount of force and displacement. However, the air pressure has to be supplied by an air pump which is rigid and the robot's body has to be tethered to be attached to the pump, which not only makes the entire robot's body ''not soft'' but also occupies large space if the source of power is taken into account.
Another type of soft actuators are based on Dielectric materials [10]- [13]. This type of elastomers are smart material systems that produce large strains. They belong to the group of electroactive polymers. Dielectric actuators transform electric energy into mechanical work. They are lightweight and have a high elastic energy density. However, they usually need extensively high amount of voltage (in order of 20KV) to be activated.
Shape memory actuators are another type of soft actuators that have been developed by researchers [14], [15]. A shapememory alloy (SMA) is an alloy that can be deformed when cold but returns to its pre-deformed (''remembered'') shape when heated. Shape memory alloys are popular as actuators for soft bioinspired robots because they are naturally compliant, have high work density, and can be operated using miniaturized on-board electronics for power and control. However, SMA actuators typically exhibit limited bandwidth due to the long duration of time required for the alloy to cool down and return to its natural shape and compliance following electrical actuation.
Liquid-gas transiting actuators [16], [17] are another type of soft actuators, where a liquid with low boiling point (usually around room temperature) such as ethanol, encapsulated inside a soft and stretchable body such as PDMS. By increasing the temperature, through generating heat (electric current through a resistive tiny wire) the liquid would change its stage to gas and as a result huge changes in volume will happen. This would lead to inflation of the soft PDMS and thus motion and force can be generated. However, this type of actuators have very limited bandwidth as it takes time to reverse the motion by cooling down the gas into liquid.
Recently, some soft actuators have been developed based on the electromagnetic working principle of traditional rigid actuators [18]- [20]. In this type of actuators, the force is generated based on the Lorentz law: a conductive element carrying electrical current in a magnetic field will experience a force acting perpendicular to both the magnetic field and direction of electric current. This type of soft actuators do not need high amount of input voltage as Dielectric actuators do, and are fast in response and can develop considerably high amount of deflection. However, the challenge in these type of actuators is the relatively low amount of generated force. The reason is due to several factors: first of all, these actuators use a conductive liquid (usually Eutectic Gallium Indium or EGaIn) which have higher electrical resistance as compared to copper wires in traditional rigid actuators. Also, these actuators usually use flexible permanent magnets made by a mixture of PDMS and magnetic particles, which have lower magnetic strength as compared to rigid permanent magnets in traditional rigid actuators. Another challenge is due to limitations on fabrication process of the micro-channels inside PDMS. Some 3D printers are able to create very small size micro-channels, however, these devices are very expensive. A Comparative analysis of fabrication methods for achieving micro-channels in PDMS is presented by [21].
We have developed a novel Electromagnetic Soft Actuator (ESA) based on the working principle of solenoids with permanent magnets [22]- [24]. In our ESA, electromagnetic field is created by applying electric current through a soft conductive coil. The coil is made of PDMS micro-pipe with diameter of around 0.1mm filled by EGaIN. In our design, two coils are antagonistically embedded inside a PDMS body with a springy connection in between as it is shown in Figure (1). Once electric current is being supplied the two coil can get magnetized based on the Lorentz law and attract each other. In order to intensify the electromagnetic field between the two coils, a flexible permanent magnet is placed inside the coils. Once electric current flows inside the channels it creates electromagnetic field. Two opposite magnetic poles attract each other and create force that deflects the spring linkage between the two coils.
The permanent magnet is made by pouring a mixture of PDMS and magnetic particles into a 3D printed cylindricalshaped mold, where the mixture in placed inside a strong external magnetic field during the curing process, to align the magnetic orientation of each particle. As a result, once the mixture is completely cured, the magnetic particles remain aligned even when the external magnetic field is removed. The result is a cylindrical flexible material with magnetic properties.
This type of soft actuator is highly scalable where the scaling factor is determined by the available manufacturing technology. For example, the smallest PDMS micro-pipe available has a diameter around 0.1mm, whereas using highend 3D printers, PDMS can be printed as a coil-shape with embedded helical micro-channel inside with diameter in the order of micrometers. However, due to extremely high price of such 3D printers, this manufacturing technology is very expensive and only justifiable in mass production.
By decreasing the size of the actuator, the generated force will be reduced. Interestingly we found that by scaling down the ESA's size, its force to volume ratio increases. This was confirmed through experiments on different sizes of ESAs [24]. This behaviour was also analytically proven using the geometry of the actuator and electromagnetic equations based on the Lorentz Law [25].
This property of ESAs suggests that by reducing the size of ESAs and attaching them in series and parallel as a network, the output force can be enhanced compared to a single ESA with the same size of the networked ESAs. Due to their light weight, their fast response, and being able to be operated with voltage range around 10V 80V, with increased output force, ESAs can reproduce a soft actuation technology suitable for pHRI, especially for rehabilitation and support of patients with mobility impairments.
As it will be shown in section 3, the equation governing the force is nonlinear and contains design variables that are both real and categorical in nature. In order to maximize the output force, therefore, we are dealing with a nonlinear mixedinteger optimization problem. We propose to use a Branch and Bound algorithm based on interval analysis technique that is capable of obtaining the global optimal with mixedinteger variables.
In order to examine the capability of a networked ESAs to be used as drive train for a rehabilitation or force augmentation device, we consider a case of an active elbow brace and use of our proposed Branch and Bound optimization framework to determine whether or not our optimal network of ESAs can achieve the required amount of force and deflection.
The organization of this paper is as follow: in section 2, the design parameters will be explained and the output force of a single actuator as well as that of a network will be formulated as a function of design parameters. In section 3, the proposed Branch and Bound design optimization algorithm based on interval analysis will be presented. Section 4, presents an active elbow brace as a case study and applies the proposed optimization algorithm to determine the capability of networked ESAs as a drive train for such a device. Finally section 5 presents the conclusion and future works.

II. DESIGN PARAMETERS OF NETWORKED ELECTROMAGNETIC SOFT ACTUATORS AND FORCE FORMULATION
In this section, we first formulate the force as a function of design parameters of a single and then double (antagonistic) coils. We then include the effect of magnetic core to calculate the resultant force. In the next step, the generated force of networked ESAs will be formulated as a function of its design parameters.

A. ELECTROMAGNETIC FORCE OF TWO ANTAGONISTIC COILS
As mentioned before, our proposed ESA consists of two electromagnetically inductive coils that contain a flexible permanent magnet core. The force generated by the coils depends on several parameters as following: number of turns, number of loops, length and diameter of the coils, electrical current, length and diameter of the flexible permanent magnet, and the type and ratio of mixture of materials being used.
To design the ESA, a theoretical analysis should be conducted on the electromagnetic force that is being produced by the two coils, taking into account that the two coils share a common flexible permanent magnetic core. For this purpose, we first calculate the magnetic field by the Biot-Savart Law [26] along the axis z of a representative loop of the coil with unit vector ofk, considering an arbitrary point P through which a steady current I is passing ( Figure (2-a)). The magnetic field at point P due to the contribution of the current element can be presented by applying Biot-Savart Law, as follow: where µ 0 is the permeability constant synonymous to the permeability of free space or as magnetic constant, − → r e is the effective radius defined as: and R is the magnitude of the coil radius vector − → r as it is shown is Figure (2-b). To calculate the magnetic field at point P, d − → B has to be integrated over the entire circular loop. The first integral vanishes due to the fact that radial unit vectors around the circle sum to zero. Hence, the overall would only be the axial component of the magnetic field as: This is the magnetic field of a single loop acting on point P along the axisk. A solenoid consists of many loops and layers of loops. In order to calculate the magnetic field acting on an arbitrary point P along the axisk in a solenoid, the resultant magnetic field of a single loop should be integrated over the entire length of the Solenoid. dz , as shown in Figure 2, is a selected current element on the solenoid. The z dimension of the selected point P is always measured from the central loop of the solenoid, where the amount of current dl passing through the element dz is nIdz , where n = N l is the coil's loop density.
The infinitesimal magnetic field due to the presence of an infinitesimal electric current within the length of dz is: Similarly by integrating d − → B z over the entire length of the solenoid, the total magnetic field at point P can be found as: This is the magnetic field of a solenoid at any given point P along its central axisk with a distance z from its center that is located along the axisk at the middle of solenoid's length.
In our ESA's design, two antagonistic coils are interacting with each other and therefore, contribution of both coils on a given point should be taken into account. By calculating magnetic field due to each coil (left and right) and then applying the superposition principle, we can find the total magnetic field at any given point P located on the common central axis of the two antagonistic coils.
To calculate each magnetic field, the distance of point P from the center of its corresponding coil (z 1 and z 2 ) should be calculated.

B. EFFECT OF PERMANENT MAGNET CORE
In order to determine the effect of flexible permanent magnet core on the resultant force, we will establish the Electromagnetic Charge Model. The Charge Model is a useful method for analyzing the force of permanent magnet. Given ρ m = −∇.
− → M (A/m 2 ) as equivalent volume charge density and σ m = − → M .n A/m as equivalent surface charge density (n is the unit vector of the surface), the resultant force due to presence of a permanent magnet can be calculated as: The permanent magnet has a consistent and uniformly distributed magnetic field along its central axis as − → M = Mẑ which makes ρ m equal to zero. The net magnetic moment per unit volume of the permanent magnet is M that can be calculated as: − → M = B r /µ 0 , where B r is remanence or residual flux density that is obtainable from the hysteresis curve of the material.
The unit vector of the surfacen has three distinct directions as follow: when z = 0,n is opposite ofẑ, when z = h (h is length of the cylinder),n is equal toẑ and when the distance from the central axis r is equal to radius of cylinder r,n is perpendicular toẑ. Therefore, if the magnetization of the magnetic core along z axis is M s , the surface charge density σ m would be −M s when z = 0 and M s when z = h.
However, of the surface of the permanent magnet where r = r, the surface charge density would be zero as the two vectorsẑ andr are orthogonal.
In our ESA design, as it is shown in Figure 3, the North pole of the permanent magnet is located in the right coil while the South pole is in the left coil. The dimension d is the distance between each coil's mid point and the magnetic poles. Each pole has the same distance to the correlated coil's end as the magnetic core is placed precisely at the middle of the line that is connecting the two coils. Due to symmetrically antagonistic geometry of our design, only the applied force at one end surface of the magnetic core has to be calculated and then doubled to obtain the total applied force. In order to find the magnetic field at any arbitrary point within this design, let's consider z 1 and z 2 to be d and d + h, respectively. The total magnetic field due to left and right side coils can then be calculated using Eqs. 4 and 5 as: where The resultant force can the be formulated using Eq. 6 as: By replacing Eq. 7 into Eq. 8, the force that applied to either North or South Pole of the magnetic core and generated by magnetic fields of both the right and the left coils can be calculated as: Therefore, the total force applied to the magnetic core both at the North and South poles would be F total = 2F. The current in this equation can be determined by the input voltage V in and equivalent resistance of the coil's circuit R eqv as I = V in R eqv . The equivalent resistance of the coil is a function of length L w and cross-section area A w of its conductive wire VOLUME 8, 2020 as well as specific resistivity of the conductive material ρ w as The total force a nonlinear function in which we are dealing with both real variables such as dimensions as well as categorical variables such as the residual flux density or specific resistivity of the materials used in ESA. Figure 4 shows components of a representative single ESA and its flexible magnetic core.

C. NETWORK OF ELECTROMAGNETIC SOFT ACTUATORS
A network of ESAs consists of ESAs attached in series and parallel, similar to the arrangement of Actin and Myosin filaments in skeletal muscles. When ESAs are attached in series fashion, they compose a Fiber and when fibers are attached in parallel, they make a Fascicle. Output force of a fiber is equivalent to the output force of a single ESA while the total deflection of a fiber is summation of deflections of its ESAs.
The output force of a fascicle, however, is the summation of output force of its fibers, which, as mentioned before, is equivalent to the summation of output force of single ESAs attached in parallel. Therefore, to maximize the output force of a networked ESAs, in addition to maximizing the output force of a single ESA, the parallel arrangement of ESAs should be optimized. This parallel arrangement is a function of the cross-section profile of ESAs. Usually coils have circular cross-section profile, but the conductive wire can also be wrapped in different cross-section shapes such as triangular, square, etc. basically if n represents number of sides of a shape, the cross-section profile can have n sides where n can change from 3 (i.e. triangle) to ∞ (i.e. circle).
The purpose of this subsection is to geometrically design a fascicle network of ESAs so that the magnetic field and force out of a certain occupied volume of actuator network can be maximized. In the previous subsection we formulated the magnetic field and output force of a single ESA mainly by taking into account the geometry of the coils embedded in it. Now, we intend to discuss the formulation of a network of soft actuators, which are positioned with their long axis parallel next to each other. Given the point that the coil is the heart of the proposed ESA, we need to study the coil's cross-section profiles. The idea is to eventually find the best cross-section profile for the axially packed network of solenoids shown in Figure 5 which yields to the maximum generated magnetic field and force out of a certain occupied volume of the network. Since, we intend to compare cross section profiles for each solenoid profile, we just need to calculate the magnetic field and force of one loop coil. The polygons and circle shapes for the coil's section profiles will be studied here. For this purpose, similar to the approach was taken in previous subsection, the magnetic field of the considered profile needs to be determined using Biot-Savart Law. For all the polygonal profiles, firstly, we need to calculate the magnetic field generated by a finite straight line representing the side of the polygons at a desired point as shown in Figure 6  The current carrying element in x direction is d − → s = dx î . the goal is to find the magnetic field of the current carrying wire at point P. The position vector describing P is − → r P = aĵ. The relative position vector which points from the source point to the field point is − → r = − → r P − − → r . for point P with the location of (0, a), the relative position vector would be − → r = aĵ − x î . The magnitude of this vector defines the distance between the current element and the point P and is equal to r = (a 2 + (x ) 2 ) 0.5 . The corresponding unit vector is given byr = − → r r = sin θĵ − cos θî. According to the Biot-Savart Law the contribution to the magnetic field due to the current carrying element is: Substituting a = r sin θ and dx = −a csc 2 θdθ in Eq. 10, the differential contribution on the magnetic field is obtained as: Integrating this infinitesimal magnetic field over the range of angles from −θ 1 to θ 2 , the total magnetic field can be calculated as: In case of symmetrical arrangement of the conductive wire around perpendicular line passing through point P (i.e. θ 1 = θ 2 = θ), the total magnetic field would be: Now that we calculated the magnetic field of a current carrying finite wire, it is possible to calculate the magnetic field produced by polygon shaped current caring loops. For each polygonal shapes, the magnetic field at the center of the section is obtainable by superposition of the magnetic field produced by each sides at the center of the profile. For instance, for a triangular profile, depicted in Figure 7 the magnetic field at the center of the section would be: where a tri is the apothem of the triangle. Therefore, in order to compare various cross-section geometries, we would equalize different sections in terms of cross-section area and find the correlated edge length of each polygon solenoid section.
In other words, we assume that all polygon sections have same cross-section area and then the side's length for each one can be calculated. The cross-section area can be calculated by the apothem of each polygon and number of the sides on each polygon as: A cross−section = m.a 2 .tan(180/m). Once, the magnetic field of all cross-section profiles has been calculated, the generated force resulted from the interaction between coil and magnetic core can be obtained: For each coil profile geometry, the magnetic core has the same profile with the correlated coil. Since, we assumed earlier that all profiles have the same cross section area, the magnetic core cross-section A cross−section would be the same as well. Therefore, the produced force from the coil would be proportional to its magnetic field for each case, hence the ratio of generated forces would be the same as the ratio of magnetic fields.
Another parameter to be considered is packing density λ that can be defined as a function of total cross-section area of the network A network and cross section area of a single ESA and number of parallel fibers within a network P as: λ = P.A cross−section /A network . If all cross-section area of the network can be occupied by ESAs, then λ = 1.
Let's say number of series ESAs within a fiber is S and number of parallel fibers within a network is P. Then the total elongation of the network would be L network = S i=1 d = Sd. The overall force that the network can exert would be: d, a, n, m, B r , I , S) (16) which is a function of the following design parameters: number of fibers in parallel, length of a single ESA, distance between two coils of a single ESA, apothem of the polygonal profile of ESA, number of polygon's side, coil density, residual flex density of the permanent magnet, the electric current and number of series ESAs within a fiber. All these parameters are real continuous variables except for number of fibers in parallel, number of polygon's side, number of ESAs in a fiber and coil density that are discrete as well as residual flex density of the permanent magnet which is categorical variables. The goal is to maximize this nonlinear mixed-integer function as an objective function subjects to functional constraints regarding the overall size (length, width and thickness) that the network can occupy, the overall deflection of the network as well as its power consumption. These functional constraints depend on particular application of the network. For this work, the functionality of the networked ESAs will be examine as actuators for an active elbow brace.

III. OPTIMAL DESIGN BASED ON BRANCH AND BOUND ALGORITHM
In the previous section, the output force of a single ESA as well as that of a network was formulated as a function of design parameters. Obviously, the force function is nonlinear and as the design parameters are both real and categorical, the optimization process should be able to deal with nonlinear mixed-integer functions. Here we propose Branch and Bound algorithm to find the optimal solution to this problem.
Generally, the problem of designing machines (e.g. actuators) is understood and formulated as an inverse problem. The direct problem of the design can be defined as follows: From an actuator where the structure, the dimensions, and the materials are known, compute the performances; for example, the velocity, the torque, etc.
In the past years, the corresponding direct problem has first been studied by developing the well-known methods such as finite-element. Thus, with these numerical tools, the designer can iteratively generate and study some possible solutions of the considered inverse problem. Nevertheless, such a way cannot be entirely satisfactory because it does not solve the inverse problem rationally. VOLUME 8, 2020 The corresponding inverse problem of design is: From the values given by the desired performance (for example the torque), get the structure, the dimensions, and the material of the required actuator. Such a problem is ill-posed in the Hadamart sense [27], because: 1) the existence and the uniqueness of the solution cannot be guaranteed and 2) this problem may generate an infinite number of solutions. Many counter-examples can be found, like the actuator considered in [28]. The inverse problem of design of electrical machines is explained in [29].
More recently, optimal dimensioning problems (that only deal with continuous variables) have been considered and solved by using classical local search algorithms, such as augmented Lagrangian method or Sequential Quadratic Programming (SQP) method, associated with analytical models [30], [31]. Many analytical models for different actuators have then been introduced depending on different assumptions [32]- [35]. Generally, in all these cases, only classical local search methods are considered and, therefore, only optimal dimensioning problems have been considered; the hypotheses used in this kind of algorithms can only deal with continuous optimization problems. Unfortunately, even for this kind of problem, it has been proven in [29] that such an inverse problem must be solved by using an exact global optimization algorithm and then it must be treated as a constrained-continuous global optimization problem.
In the classical literature about modeling and optimal design of actuators, some particular inverse problems, named predimensioning problems of design, were considered [36]. These works focus on the association of analytical models with local standard optimization algorithms [37]. Thus, first optimal solutions are found for some actuators. In this case, the solutions obtained by local optimization methods depend on the starting point introduced by the designer. Stochastic global optimization algorithms are not well adapted to solve this kind of problem because there are too many hard constraints. In [30], it is proven that even if a simple actuator is considered, local search solutions failed to find the optimal solution; see [36] to compare results with the exact ones published in [30].
In [32], [38], a rational way associating combinatorial analytical models of actuators and an exact global optimization algorithm, was proposed and studied. The algorithm was an exact global optimization method; see [30], [33], [35] for details on such algorithms and on their rigorous convergence to the global optimum. The obtained solutions satisfied the imposed set of conditions via combinatorial analytical models. When these optimal solutions were validated by the means of numerical tools (such as finite-element methods), some differences about the values of the electromagnetic torque were noted. Thus, the optimal results are not robust and this involves some adjustments of the parameters of the obtained solution. These adjustments can be done by iteratively solving the direct problem of design until the desired performance will be satisfied using a numerical model. In some numerical tools, such as ANSYS, stochastic and local optimization algorithms can be used to solve this problem of optimal adjustments.
The principle of this algorithms is based on subdivisions of the considered initial domain into smaller and smaller parts such that one can determine, using interval analysis [39], which box can be discarded. This is because interval computations will produce the proof that a box cannot contain the global solution or prove the fact that at least one constraint could not be satisfied in a box. Therefore, at the end of the algorithm, the global optimum will be enclosed with a given accuracy, however the algorithms must be capable to deal with both real continuous and categorical variables.
The main idea is to subdivide the initial domain space X ⊆ IR n × c n i=1 K i into smaller sub-boxes Z ⊆ X and to delete the considered box Z , if and only if it can be proven that Z cannot contain the global optimum.
In order to subdivide the domain into smaller boxes, we need to consider the space in which the boxes are located in. In this problem, this space is a combination of continuous as well as categorical variables.
Branch and Bound algorithms have been modified to deal with integer variables, such as in our previous work in finding an optimal solution to the facility layout problem [40]. The categorical variables cannot be directly considered in subdivision of a box and an expression of a function. Here, we propose a method to change categorical variables into integer ones. This is done through an univariate function a that assigns an integer to a categorical value. For example, if one material has flux density of 0.85 Tesla and stiffness of another material's flux density is 0.9 Tesla, this univariate function assigns integer value = 1 to the first stiffness value and integer value = 2 to the second one: a(1) = 0.85 Tesla and a(2) = 0.9 Tesla. Therefore, while integer values 1 and 2 can be used for the purpose of subdividing the boxes, the values a(1) and a(2) will be used in calculating the objective function and constraints.
The classical principle of subdividing is to choose a coordinate direction parallel to which Z has an edge of maximum length. Then, Z is subdivided normal to this direction [39]. For continuous variables the subdivision of a box Z on its kth component into to boxes Z 1k and Z 2k will be applied at the midpoint (as defined in [39]) of the original box Z as follow: , where Z L k and Z U k denote the lower and upper bounds of kth component of box Z , respectively. For categorical variables, however, this subdivision rule has to be slightly modified as follow: where [x] I is the integer part of x. At each step of the Branch and Bound algorithm, and as the branches are generated through subdividing the boxes, we will have an interval for each variable. The interval analysis proposed by [39] is a powerful tool to calculate upper and lower bounds of a function over a box. This is simply done by the concept of natural extension of a function.
The natural extension of an expression of f into interval consists by replacing each occurrence of a variable by its corresponding interval (which encloses it), and then by applying the rules of interval arithmetic as explained in [39]. For example, if f (x) = x 2 + x + 1 and x ∈ X = [1, 2], then F(X ) = ([1, 2]) 2 + [1, 2] + 1 (X and F(X ) are both intervals) is the natural extension of f (x) over box X . It has been proven by [37] that the natural extension is always an inclusion function, which means that it encloses the upper and lower bounds of f over the box (or even all boxes when dealing with multiple variables) X . Therefore, Interval arithmetic is only defined for continuous functions, and thus in our method, the inclusion functions must be extended to deal with discrete variables as well. As mentioned before, we will convert categorical variables into integer variables through univariate functions. These integer variables will then be further relaxed into continuous variables. For example, if an integer variable belongs to set Z L , Z L + 1, Z L + 2, . . . , Z U , a continuous interval [Z L , Z U ] will be considered for this variable. Replacing an integer set by its corresponding relaxed continuous interval, an inclusion function can then be constructed. The proof is obvious because it comes from the fact that the relaxed compact interval sets enclose by definition the initial discrete sets.
In order to proceed with the elimination of the boxes that cannot contain the global optimum solution, just the computation of lower bound of a given function f over a box Z , i.e. lb(f , Z ), is needed. We will compute these bounds by using interval analysis as follow: Considering thatf denotes the current solution (in fact, it is just the best evaluation of f at this stage of the algorithm such that all the constraints are satisfied), one obtains the following: 1) No global solution is in Z , if ib(f , Z ) >f , a lower bound of f over Z is greater than a solution already found, then no point in Z can be a global minimum.
2) No feasible solution is in Z , if it exists k such that ib(g k , Z ) > 0, or it exists k such that ib(h k , Z ) > 0 or ub(h k , Z ) > 0 (ub = upper bound). In this case, a very small positive real value can be considered instead of 0, in order to address numerical approximation difficulties due to the floating point operations [39].
This implies that at the end of the algorithm, accurate enclosures of the global minimum value and of all their corresponding solutions can be expected. Indeed, at each step of the algorithm, one has the following properties for all the remaining sub-boxes Z ⊆ X : 1) lb(f , Z ) ≤f and 2) all constraints are always satisfied.

IV. CASE STUDY: AN ACTIVE ELBOW BRACE
As mentioned before, bio-inspired arrangement of ESAs can lead to an enhanced output force. If we can achieve the level of forces that are exerted physiologically at a human's extremity joints, by on-board batteries, this actuation technology can greatly impact rehabilitation or force Algorithm 1 Branch and Bound Algorithm 1: set initial box X := IR n × c n i=1 K i 2: set current solutionf = +∞ 3: set initial lower bound = (+∞, X) 4: extract from the lowest lower bound lb 5: divide the chosen box V by its midpoint m to get subboxes V 1 and V 2 6: for j = 1,2 do calculate v j = lb(f , V j ) calculate lower bounds of all constraints over V j using interval analysis iff ≥ v j and all constraints are satisfied then: if and only if all constraints are satisfied at m j c: iff has been changed from its previous value, then remove all (z,Z) from where z >f 7: if min (z,Z )∈ z =f then STOP, else GoTo step 4 8: end if 9: end for augmentation applications for mobility impaired patients, as it is light weight, has high bandwidth, and it is powerful and portable.
To determine whether or not networked ESAs has this capability, we selected a human elbow joint, as a critically important joint in performing our daily activities. The elbow joint functions as a functional link between the upper arm and the forearm. It positions the hand in space, allows the forearm to act as a lever in lifting and carrying, and provides precision in both open and closed kinetic chain work. Therefore, even a mild impairment in elbow can significantly reduce the ability of the hand to reach its objectives, with the problem being compounded by the fact that the other joints in the upper limb are unable to compensate for this loss.

A. CONSTRAINTS
In this section, we first determine the performance requirement of an active elbow brace that is supposed to be worn by a patient and match the performance of a healthy elbow joint [41], in terms of maximum torque and angular deflection, i.e. flexion-extension. Range of elbow's flexionextension in a healthy individual is between 0 to 104 degree, while the maximum torque, that the elbow can handle at 90 degree (flexed) is about 4-6 Nm [42]. The other performance parameters of the active elbow brace that are summarized in Table 1, have been calculated based on available passive elbow braces in the market.
We consider a popular elbow brace by Orthomen [43]. To make this passive brace active, two networked ESAs can be attached to the brace in an agnostic-antagonistic fashion, similar to the Biceps and Triceps skeletal muscles. Hence, when one network deflects, the other one relaxes to perform flexion, and vice versa for extension. The antagonistic ESA network is split into two parts with a rubber in between that can slide on a roller attached to the brace. Each ESA network is pinned to the brace structure so that it allows for rotation. With this setup, the agnostic ESA network will have a length of 0.24m when it is relaxed and 0.15m when it is fully deflected. That means the total elongation of the network would be 0.09m. Therefore, the overall length of the network L netwrok = S i=1 l ESA = S.l ESA should be less than 0.24m. Based on structure the brace, the cross-section area that the network can occupy is estimated to be less than 0.24m and L network equal to 0.09m. Furthermore, considering the available space around the brace, the overall thickness T network and width W network of the network are up to 0.02m and 0.05m, respectively. This would lead to a maximum allowable volume of the network V network around 10cm 2 .
At elbow angle = 90 degree, the level arm (distance from force vector to the elbow's center of rotation) is about 0.0845m. Considering 40Nm torque, that means the overall force from the agnostic network should be about 71N. Figure 8 shows a conceptual schematic of such a design. The optimal force of the network will be compared to this value.

B. RANGE OF PARAMETERS
In addition to the constraints, ranges of the input parameters have also to be determined such as how much electric power is available to the network in terms of Voltage and Current. Since the active brace is supposed to be activated using an on-board battery, we set the maximum available power as 10 Watts (average output of a power bank). Length of a single ESA should be limited to the maximum length that a network can have, i.e. 0.24m. This is the case when each fiber is composed of only one ESA. Similarly, the distance between the two coils should be limited to 0.09m. The apothem of the cross-section area should be limited to 0.02m which is the thickness limit of the network. We limit the maximum number of sides of the cross-section polygon to 8, however to consider circular cross-section profile, we would also assign a large number (e.g. 100) as for the polygon's number of sides. Coil loop density is limited to 5000 loop per meter, due to the smallest feasible diameter size of the flexible PDMS wire (i.e. 0.1mm). The categorical variable, i.e. flux charge density has to be experimentally quantified. As mentioned before, to prepare the flexible permanent magnet, a mixture of PDMS and magnetic particles were cured inside a 3D printed mold while the whole mixture was exposed to a strong external magnetic field. However, the ratio between magnetic particles and PDMS plays a very important interdisciplinary role: i.e. both in magnetic and mechanical domain. By increasing the ratio of magnetic particles, the flux charge density increases that would lead to enhancement of the magnetic field and the force generated by the coils. However, in mechanical domain, this increase in the magnetic particle ratio, affects the Young Modulus (i.e. elasticity) of the flexible magnet as well as its yield point, which in turn, affects the maximum passive elongation ratio of the permanent magnet. We prepared 5 flexible permanent magnet samples each with different mixing ratios (between 8% to 28% with incremental mixing percentage of 4% of magnetic particles). For each sample, the flux charge density was measured by Magnetic Field Instrument (MFI), a device used to measure the magnetic field or flux around permanent magnets, coils, and electrical devices. The flux density for each mixing ratio was presented in an interval domain as the mixing ratio could not be precisely adjusted. The result is shown in Table 2. The optimization problem then can be formulated as: n ∈ 1, .., 5000 P ∈ 1, .., 100 S ∈ 1, .., 100 B r ∈ 1, .., 5 : Table2 m ∈ 1, .., 8and100 (17) where P r ∈ IR 4 and P c ∈ 5 i=1 are real and discrete variable vectors, respectively. With these settings, we performed the proposed Branch and Bound optimization process to find what design parameters would lead to maximum output force and determined if that optimal design parameters can satisfy the aforementioned performance and input parameters.

C. OPTIMIZATION RESULTS
We conducted mechanical tensile tests on several mixtures of PDMS and magnetic particles with different mixing ratios. The results suggested that mixing ratio of magnetic up to 28% (with respect to the weight compared to that of PDMS) would lead to safe (i.e. passive elongation ratio > 35%).
Result of the optimization problem is presented in Table 3. The Results on the optimal cross-section profile showed that although for a stacked network of ESAs, the triangular cross-section profile makes the largest magnetic field and force at the center of coil compared to the other profiles, it also consumes the most electrical power such that after normalizing and obtaining the ratio, the triangular profile is the least efficient section profile while the circular one is the most.
We also considered the effect of packing density of the profiles since all of the coil's section profiles do not cover the available space entirely and there are some free spaces among the profiles for circular, octagonal and pentagonal arrangements. For this purpose the packing density coefficient (λ = P.A cross−section /A network ) was defined and calculated for each section profiles. After applying the packing density coefficient to the previous findings we found that the hexagonal arrangement for the coil's section profile (m = 6) is the most efficient one due to have the largest amount (λ = 1) among all section profiles.
Furthermore, the optimal agnostic network on each side (left and right) should consist of 40 fibers, attached in parallel, each has of 9 ESAs in series for the agnostic networks. With these networks the total weight of the active brace would be 200 g, while the volume stays below 200 cm 3 and the flexion-extension range achieves 93 degrees rotation around the elbow joint (11 degrees less that what is actually required). The force of the agnostic network as a result of optimization is 71 N which means that the output torque of the active brace would be 6 Nm that matches the required output torque at the elbow according to [42].
In order to evaluate the result of the optimization process, we have built a small network based on the optimal design parameters with 3 parallel fibers where each fiber consists of 4 series ESAs as it is shown in 9. The small network was able to lift a 2.5N weight, with its 3 fibers. When scaled up to the number of fibers of agnostic network suggested by the optimization process, the total output force would be 66.6N, which is slightly less than what was predicted by the optimization process (i.e. 71N). A reason could be due to actual coil density. As it is clear from Fig. 9, there are gaps between the micro-channels inside an ESA, as the flexible wire was manually wrapped around the coil. This would lead to less number of the loops per meter and eventually less output force. An automated wrapping technique would solve this issue.

V. CONCLUSION AND FUTURE WORKS
In this paper, optimal design for a network of novel ESA based on branch and bound algorithm was presented. The novel electromagnetic soft actuator operates based on the working principle of Solenoids, that consists of two antagonistically located coils made of flexible wires and share a flexible permanent magnetic core.
It was shown that by reducing the size of these ESAs, the force to volume size ratio increases, which suggest a network of miniaturized ESAs would achieve higher amount of output force compare to a single ESA with the same size of the whole network. In this work this was numerically tested for a case study of an active elbow brace.
Branch and Bound optimization algorithm was employed to achieve optima design of a single ESA and consequently optimal design of a networked ESAs to achieve the maximum output torque at the elbow joint, while the performance parameters are being satisfied.
The result showed that having a network of ESAs as drive train for an active brace, we can satisfy the performance parameters, for supporting the elbow joint of a patient with decreased muscle performance and mobility.
This suggests that with the available manufacturing process discussed in this paper, the actuation technology based on electromagnetic soft actuators can to be used as drive trains in robotic prosthesis and robotic exoskeletons, to support patients with decreased muscle function at their affected joints. Our future endeavors are focused in enhancing further the produced torque by the ESAs, in order to be utilized for robotic prosthetics and exoskeletons in patients with complete loss of muscle function.
This actuation technology is uniquely suitable in rehabilitation and/or force augmentation applications for those mobility impaired patients that have not completely lost the ability to move their affected joints and would need some extra help to recover or be able to perform their daily tasks. Considering the huge population of these types of mobility impaired patients (e.g stroke patients, peripheral arterial disease, traumatic injuries, neuropathies, senescence and frailty) electromagnetic soft actuators provide novel potential solutions, for wearable and next-to-skin type of assistive technologies, at low production cost, safe, portable and yet sufficiently powerful with low power requirement and high bandwidth.
Planning for the future, we will manufacture the active brace powered with the network of ESAs as suggested by the Branch and Bound algorithm and test its potentials in rehabilitation trials, such as those for elbow stiffness or isokinetic motion for elbow spasticity and other motor dysfunctions from various medical pathologies.