1. Introduction
The energy consumption of industrial machinery is a topic of primary importance due to environmental and economic considerations [
1]. The 45% share of electric motors in the global electric consumption [
2] supports the statement that any energy-saving method should be investigated thoroughly. The methodology proposed within this paper applies to all planar four-bar mechanisms with an imposed movement of the end-effector and/or output link BC (see
Figure 1). The potential of this scope was indicated in [
3,
4,
5], stating that four-bar linkages are extensively used in practical engineering applications. Moreover, reciprocating Point-To-Point (PTP) machinery is progressively common within the industry [
5].
A required PTP displacement of the output link BC can be obtained with different link lengths. Therefore, the geometry parameters depicted in
Figure 1 can be considered as design parameters (DPs) in a mechanism’s optimization. The design optimization of a PTP mechanism can reduce the energy consumption (linked to the RMS motor torque) of electric machinery, as indicated in
Figure 2. Awareness about the influence of machine components’ geometry on energy consumption has recently attracted attention [
6,
7,
8]. Mechanism models [
9,
10] replace prototyping, allowing the computational evaluation of multiple designs with limited costs.
An emergency ventilator was used as a validation case within this study. This mechanism was constructed during the first wave of the COVID-19 pandemic by a non-profit organization [
11]. Having continuous (24/7) access to electricity is not obvious within low- and middle-income countries. Thus, having a mechanism that consumes a minimum of electric energy enabling the usage of batteries is highly relevant. Therefore, the objective of this study was to find the optimal design (lengths
,
and
) leading to a minimal
for a reciprocal four-bar mechanism.
State-of-the-art techniques generally use heuristic optimizers that cannot guarantee finding the global optimum [
12]. However, the method introduced in this paper guarantees the global optimum in a build’s convex hull by combining CAD motion simulations with sparse interpolation techniques. Once the constraints that limit the design space were defined, CAD motion simulations were used to sample the objective function. As such software is well known in the machine designer community, this approach allows broad industrial applicability by avoiding the cumbersome analytic derivation of the mechanism’s dynamics. Furthermore, sparse interpolation avoids the infeasible computational burden of numerous CAD simulations. This interpolation method developed by the authors and introduced in this paper is described in depth. Moreover, a step-by-step procedure is presented to allow its application to any mechanism without the necessary in-depth insights into the underlying mathematics. Advances in sparse interpolation [
13] make this approach very scalable to higher-dimensional optimizations. The interpolation model delivers a mathematical description that enables global optimizers such as INTLAB [
14] even in high-dimensional problems. CAD motion simulations avoid cumbersome modeling, making this method also scalable to more complex mechanisms. In conclusion, one can state that the scalability of the proposed approach is enabled by the novel interpolation technique, which limits the number of CAD motion simulations to an absolute minimum.
In the literature [
15], minimizing the driving torque of simplified mechanisms is accomplished by utilizing dynamic equations to predict the system dynamics. However, as noted by [
5], the analytical computation of mechanism dynamics can become complex when dealing with systems containing a large number of components, as is commonly the case in industrial machinery. Alternative approaches, such as the use of generic equations of motion, have been proposed in the literature [
16]. However, these methods still typically require information such as the reduced motor inertia as a function of motor position, which can be obtained by the method of kinetic energy [
17], being cumbersome and error-prone. More specifically, studies [
18,
19] on generic dynamic equations for four-bar mechanisms still necessitate the determination of the center of gravity position and mass for each component of the mechanism with every design modification. This presents a challenge for machine builders, as it requires extensive case-specific analysis and data collection. Therefore, in engineering practice, it is standard to use Computer-Aided Design (CAD) software to build a dynamic model.
Moreover, References [
18,
20] did not define the feasible search domain nor include it in searching for the optimum result. The constraints that define the feasible design space are important as defects, giving infeasible designs [
21], frequently occur in the kinematic mechanism synthesis of a four-bar linkage. The optimization algorithms of [
15,
18,
22] ensure that the objective function converges towards a minimum, yet it is generally not guaranteed that the designed linkage will be feasible. Therefore, the necessary constraints should be added so that the optimal solution can fulfill the movement without inconveniences. A constrained global optimization algorithm requires a deterministic mathematical description of the constraints to find the global optimum. To the authors’ knowledge, this has not been done yet in the literature [
23].
Developing a four-bar mechanism that follows the desired output trajectory is a classic design problem that researchers have extensively explored [
24,
25,
26,
27,
28]. However, all the methods above are not implementable in global optimizers as the algebraic expression (when provided) is only evaluated in discrete defined points
on the coupler curve
B(
) (shown in
Figure 1). Thus, these cannot deliver a deterministic mathematical description of the feasible design space, which is required.
This paper shows how CAD-based motion simulations combined with a novel sparse interpolation technique enable a global optimizer that guarantees the global optimum strictly within the convex hull covering the larger part of the feasible design space and thereby outperforming heuristic optimizers regarding energy savings.
The mechanical design of systems is mainly performed in CAD software. These CAD models include all required information (i.e., volume, mass, friction, damping, joints, etc.) to model the dynamics of a mechanism. This information is necessary to calculate the necessary torque of the mechanism through motion simulations. The required motor torque to move the mechanism is compulsory to evaluate a certain design with all attendant external forces on the mechanism. By driving the mechanism with the motion profile
at Point O (
Figure 1), the location where the mechanism is driven in reality by a motor, the user can extract the necessary torque from the software (as in
Figure 2) to fulfill the prescribed movement
of the output link BC. Furthermore, within these motion simulations, the design parameters
,
, and
of the four-bar mechanism can be parameterized to simulate different designs. The objective value to be minimized by the optimizer is the RMS torque (
) value, necessary to drive the mechanism fulfilling an imposed PTP motion (
). The literature states that minimizing the
corresponds to reducing the energy losses in the system [
5].
Hence, by calculating the RMS torque based on CAD simulations as elucidated in
Section 2, the objective value for a certain design (i.e., certain values for the three design parameters
,
, and
) is obtained. The whole simulation process to obtain the objective value for different design parameter combinations (
,
, and
) is automated, as elucidated in [
29]. Constraints on the design parameter values are necessary to define an area containing feasible designs, as discussed in
Section 3, from which designs are selected to simulate their corresponding objective value (
). The proposed approach requires a kinematic analysis of the mechanism to find the constraints. However, this kinematic analysis is much more straightforward than modeling the necessary torque to drive the machine. Kinematic analysis is possible based on a graphical representation of the mechanism, as only the criteria described at the end of
Section 3 have to be fulfilled. Modeling the necessary torque, however, requires defining position variable inertias, friction, possible non-linear external load forces, etc.; this is a severe hurdle for machine builders and an error-prone approach. All this is solved by relying on CAD motion simulations.
Computational simulation time becomes a burden as one design evaluation can take 1 min and 25 s on average. Therefore, a wise selection of the simulated designs within the feasible design space is essential. The brute force method requires an inconceivable number of
motion simulations, with
g being the granularity of the sampling and
d the number of design parameters. Even with state-of-the-art interpolation techniques [
30], the construction of the objective function would require at least
samples, with
n the total number of terms in the mathematical description of the objective function. In the case of the emergency ventilator, this would mean 782,933 samples are required. Therefore, the selection of samples is performed with certain rules in order to use an innovative multidimensional sparse interpolation approach [
13]. This novel interpolation procedure, introduced in
Section 4, allows obtaining a model of the objective function with a sparse sampling method within the feasible design space. This reduces the number of required samples to 618, with an additional 1252 validation samples. Limiting the number of samples (CAD motion simulations) to a bare minimum to model the objective function is a major enabler for a global optimizer. In this case, the number of necessary samples is reduced from 10,000,000 to 1870. While
Section 4 discusses the interpolation technique and its underlying mathematics in depth, applying this technique only requires following the steps summarized at the end of
Section 4.
2. CAD Motion Simulations
In kinematic analysis, the linkage dimensions
,
,
, and
are known, and the resulting output motion
(and its derivatives) can be calculated. On the other hand, dimensional synthesis is regarded as the inverse, in which, for a specific output motion
, the feasible dimensions of the linkages are obtained [
31]. This paper is based on the dimensional synthesis of a planar four-bar function generation [
32]. As shown in
Figure 3, the movement
of the output linkage BC caused by
is described by a starting angle
and an end angle
. In this paper, the machine designer only defined an output motion
, which results in a reciprocal movement between the positions
and
.
The validation case is clarified to make all the following more tangible. This mechanism, shown in
Figure 4, can ventilate a patient by pressing the indenter into the bag, which causes airflow toward the patient.
Figure 4 presents the CAD model of the emergency ventilator and illustrates that the red beam, connected to the indenter (i.e., the end-effector), moves by rotating Input Link OA around Point O. This is the point where an electric motor drives the mechanism. The red beam has two predefined angles: an angle
that holds the mechanism in a position where the indenter touches the bag and an angle
that corresponds to a position in which the air is compressed out of the bag.
Figure 4 clearly shows that the mechanism is a four-bar linkage on which the method proposed in this paper can be applied.
However, the output motion of the emergency ventilator is described by the angle
(linked to the red beam), while the four-bar mechanism has an angle
that is linked to the output link BC. The relation between these angles is stated as
where
allows a conversion from
to
. The parameters
k and
in Equation (
1) (
Figure 4) are constant values that change neither in the optimization, nor during the four-bar mechanism’s movement.
A CAD motion simulation [
33] can determine the necessary torque to drive the mechanism at Point O only if the required position profile
is known at that Point O. However, the user solely defines the required position profile of the end-effector, in this case
. According to Equation (
1), we obtain
. The conversion of
to
depends on the values of the design parameters
,
, and
. Therefore, each selected design was analyzed by two motion simulations, as indicated in
Figure 5. If the design is combined with the required output motion
, the first kinematic motion simulation can extract the required motor position displacement
. Subsequently, the motor motion profile
is used in the second motion simulation, determining the required driving torque. This process with kinematic simulation and subsequent torque calculation extracts the objective value for predefined designs, which is an approach that was explained in detail in [
29].
4. Multidimensional Sparse Interpolation
Determining the objective values using CAD multi-body simulations, for the three-dimensional design problem, in a brute force way would lead to a tremendous computational burden requiring 10,000,000 samples. Therefore, we relied on a sparse data fitting method to determine a mathematical model for the objective function
. This novel interpolation technique is now explained in detail, yet be aware that conveniently applying the method can be performed using the steps summarized at the end of this section. While several interpolation methods are characterized by a trade-off between model accuracy and computing cost, sparse interpolation does not involve such a compromise. The technique introduced here uses a divide-and-conquer approach [
13,
40] by splitting up the involved numerical linear algebra problems into smaller and better-conditioned independent sub-problems.
For this method, the objective value, within the feasible design space as defined in (
7) and (
8), is determined on
l distinct lines in 3D space that are all parallel with a chosen vector
. We let
indicate the 3D vector that the
i-th parallel line is shifted over with respect to the line through the origin spanned by
for which we take
. Then, the equidistant samples on these parallel lines, as depicted in
Figure 13 (left), are denoted by:
As indicated further below, a translation of the domain to include the origin has no effect on the procedure, as this effect is at each step absorbed in the coefficients. Let us compactly denote the tuple of design variables
by
and let
denote the standard inner product in 3D space. On each
i-th parallel line, the samples
can be modeled by the sparse interpolant:
satisfying
Note that the effect or influence of is absorbed into the coefficients in , which models the behavior of on the i-th line.
The model for
can be computed using any of the existing 1D exponential fitting methods, such as [
40,
41,
42,
43]. The number of terms
in the sparse model can differ on each
i-th line. The
l individual models are only valid on their respective line spanned by
and shifted over
. Now, we need to blend these individual sparse models into an overall sparse model, valid in the convex hull of the
l lines (blue area in
Figure 13, right), which should cover the larger part of the region of interest. This requisite actually dictates the more proper choices for
and
.
In what follows, we considered every design parameter combination
U in 3D space to lie on some line parallel with the one spanned by
, also if
U is not an interpolation point. All the points on such a line take the form:
The normal plane through the origin and orthogonal to
is given by the equation:
or more compactly:
The intersection point
R of the normal plane with (
10) is thus given by
or more explicitly:
Hence, the distance of
U to this intersection point
R, expressed as a multiple of
, equals
and the points on the line given by (
10) can be re-expressed as
On each line through a point
U parallel with
, the intersection point
R with the normal plane, shown in
Figure 14, is where
on the line.
Therefore, we propose a blended 3D model of the following form to represent the overall objective value
:
where the parameters
and the value of
are already determined and where, furthermore, the overall model continues to interpolate the values
in the sample points
. The blended model (
13) coincides with the 1D models (
9) on each parallel line, and in between the lines, the exponential terms fade in and out. Since
this means
In other words, on each data line, the model consists of only terms, while in the convex hull of the parallel lines, it consists of terms. Remember that all of l and are small integer numbers.
From Equations (
9) and (
14), we consequently find
Note that remains constant along each line of the form and only varies with the projection R of that line onto the normal plane.
It remains to determine
. These functions are determined by the interpolation conditions given in Equation (
15). A simple model for
is a 2D polynomial interpolant
, as we outline now. Let us denote the intersection point coordinates given in Equation (
11) by
. The collection of points on a particular line perpendicular to the normal plane, say, here, through
R, is entirely identified by the remaining two degrees of freedom that pinpoint the intersection point of such a line with the normal plane. Since every point
on the line perpendicular to the normal plane and passing through
R satisfies the conditions:
we can take any two of the values:
to characterize the full line. Over the whole of such a perpendicular line, the right-hand sides of Equation (
16) are constant and independent of the points
U on the line. The right-hand sides of Equation (
16) are only determined by
and
R. Say, for now, that we take the first two of (
16), without any loss of generality:
. For the
l parallel lines on which the samples were collected, we find
Let us abbreviate the values in the right-hand side of Equation (
15) by
and replace
in Equation (
13) by the more appropriate
, since the interpolation conditions for
hold for a whole line and vary only with the intersection point of such a line with the normal plane:
Finally, the 2D polynomial interpolant:
where
denotes the well-known Chebyshev polynomial (of the first kind) of degree
n, is computed from the interpolation conditions:
We now apply the above to our four-bar problem. The region of interest for the design variables
and restricted by the conditions (
7) and (
8) is shown in
Figure 15, and the sampling performed in this region is shown in red in
Figure 15 (right). We took
and
to guarantee maximal coverage of the region of interest. Furthermore, the whole domain is translated over
to start sampling at the origin, in line with our description. Only 618 samples were determined by the simulations explained in
Section 2, which shape the objective function. We found that
for all
, thus yielding
terms in the global model
. The coefficients
were interpolated by a linear combination of the seven bivariate Chebyshev polynomials
and
. As a final step, we validated the blended model by collecting 1252 more simulation data on 10 other lines within the convex hull, along directions different from
. These evaluation directions are shown in purple in
Figure 16 (left), and the result of this validation is shown in
Figure 16 (right). In
Figure 16 (right), the red and purple markers depict the simulated data, and the blue markers represent the value computed from the blended model (
17). Each partial curve shows the function values of
restricted to one of the lines where the samples were collected, either for interpolation (red) or validation (purple). The overall Root-Mean-Squared Error (RMSE) equals 0.0281Nm, indicating a very good fit. When restricting our attention to
values below 5—reasonable to locate a minimum—the RMSE reduces to 0.0153 Nm. Therefore, one can conclude that the global optimum within the convex hull can be found with an RMSE for the accuracy of 0.0153 Nm, being negligible.
After this validation, we looked for a minimum of the modeled
(
17) in the convex hull of the parallel lines shown in
Figure 15 (right). This was fulfilled through a brute force evaluation of the objective function (
17), in which 10,000,000 calculations were performed in 3 min. Thus, we can conclude that the most-time-consuming was the collection of all 1870 samples, in which generating a sample took 1 min and 25 s of simulation time on average. All aforementioned timings resulted from simulations conducted on a six-core Intel Core i7-9850H with 16GB RAM.
In summary, the interpolation model within this paper was obtained through the execution of Algorithm 1.
5. Conclusions
This paper proposed an industrially applicable and scalable approach that guarantees the global optimal design, strictly within the convex hull of the sampled lines, of a four-bar mechanism based on CAD motion simulations and sparse interpolation. The process of sampling the objective value
for a combination of design parameters
,
, and
was automated by employing CAD multi-body motion simulations [
29]. Through the utilization of the motion simulations, one does not perform the error-prone derivation of the cumbersome dynamic description of any mechanical system to obtain the motor torque. Subsequently, the constraints limiting the feasible design space were introduced based on the position analysis of the four-bar mechanism. These constraints guarantee that all designs the optimizer considers are assemblable and no circuit or branch defect will occur during the mechanism’s movement.
The model reached a minimal value of 2.5989 Nm at
,
, and
. If the unconstrained design space were considered by a brute force approach, 10,000,000 objective value samples would be required to achieve the same result. As each objective value sample requires a simulation of approximately 1 min and 25 s, this would be practically impossible and seriously hamper the identification of the global optimum. However, thanks to the mathematical description of the design space constraints introduced in this paper, sparse interpolation can be applied. The innovative sparse interpolation technique described and applied here reduced the number of necessary simulations to only 618. The sparse interpolation method can be applied without an in-depth mathematical background, as it was developed as a step-by-step procedure requiring user input. This allowed us to identify the global optimal design within the obtained convex hull covering the larger part of the feasible design space. As shown in
Table 1, the method clearly outperformed the best result (local optimum) obtained through the HEEDS Sherpa heuristic optimizer [
44] and other attempts conducted in [
29] to find a better result with broadly used optimizers. The obtained global optimum was 38% more efficient than the local optimum. Moreover, this result also reduced
by 67% compared to the original design, which means that the mechanism can operate with a smaller, and thus less-expensive, motor.
The results of this study confirmed the effectiveness of sparse interpolation in designing optimal mechanisms and suggested that this approach can be applied to more complex models in the future. To this end, we intend to extend our work beyond four-bar mechanisms, leveraging the scalability and flexibility of the CAD methodology and sparse interpolation, which have no inherent limitations in terms of model complexity. Our goal is to further optimize the optimization procedure by removing the need for prior knowledge of the constraint design space, which would eliminate the need for mechanism analysis.