Faster Run-to-Run Feedforward Control of Electromechanical Switching Devices: a Sensitivity-Based Approach

Faster Run-to-Run Feedforward Control of Electromechanical Switching Devices: a Sensitivity-Based Approach* Edgar Ramirez-Laboreo, Eduardo Moya-Lasheras, and Eloy Serrano-Seco Abstract— Electromechanical switching devices, such as solenoid valves, contactors, and relays, suffer from undesirable phenomena like clicking, mechanical wear, and contact bounce. Despite that, they are still widely used in industry due to their various economic and technical advantages. This has encouraged the development of controllers aimed at reducing the collisions that occur at the end of the switching operations. One of the most successful approaches has been the use of iterative techniques. However, these algorithms typically require a large number of operations to converge, which is definitely a clear drawback. This paper presents a strategy to improve the convergence rate of such controllers. Our proposal, which is based on the sensitivity of the control law with respect to the parameters, assumes that the performance of the system is more heavily affected by some parameters than others. Thus, by avoiding movements in the directions that have less impact, the search algorithm is expected to drive the system to nearoptimal behaviors using fewer operations. Results obtained by simulation show significant improvement in the convergence rate of a state-of-the-art run-to-run feedforward controller, which demonstrates the high potential of the proposal. I. INTRODUCTION Solenoid valves are commonly used in hydraulic circuits [1] or internal combustion engines [2], while electromechanical relays can be found in power circuits, medical devices, automotive applications and, in general, in virtually all industries [3]. These are just a few examples of the widespread use of electromechanical switching devices. The open-loop dynamics of all of them is similar: when energized, a magnetic force attracts a moving component of the device, causing it to continuously increase its speed until the end of the stroke, where a violent impact occurs, generating acoustic noise and gradually wearing out the device. Despite this drawback, no alternatives have yet been found that can compete with the very low cost and numerous electrical and mechanical strengths of these devices. For this reason, research on modeling, analysis, and control of electromechanical switching devices is still being published. In order to reduce the undesirable phenomena associated with impacts, several control schemes have been proposed that aim to reach the final position with zero velocity. * This work was supported in part via grants PID2021-124137OB-I00, TED2021-130224B-I00, and CPP2021-008938, funded by MCIN/AEI/ 10.13039/501100011033, by ERDF A way of making Europe, and by the European Union NextGenerationEU/PRTR, in part by the Government of Arag´ on - EU, under grant T45 23R, in part by the “Programa Investigo” funded by the European Union - Next Generation EU, and in part by Fundaci´ on Ibercaja and the University of Zaragoza, via grant JIUZ2023-IA-07. The authors are with the Departamento de Informatica e Ingenieria de Sistemas (DIIS) and the Instituto de Investigacion en Ingenieria de Aragon (I3A), Universidad de Zaragoza, 50018 Zaragoza, Spain, {ramirlab, emoya, eserranoseco}@unizar.es This is commonly known as soft-landing control. Different approaches are reported in the literature: backstepping control [4], sliding-mode control [5], extremum-seeking adaptive control [6], or iterative learning control [7], among others. Given the fast and nonlinear dynamics of these devices, any controller can benefit from model-based feedforward terms [8], as these generally allow for improved response time and tracking accuracy. In extreme cases, where neither the position nor any other related variable can be measured for technical or economic reasons, the controller could even be based solely on the feedforward term and thus implemented in an open-loop fashion. There is however an obvious problem with this strategy: any difference between the dynamics of the real system and that of the dynamical model used to design the controller will result in a loss of performance. These differences may be due, e.g., to inconsistencies in the model equations or errors in the estimation of its parameters, or simply changes in the system dynamics due to wear or varying ambient conditions. One way to minimize the problems associated with feedforward controllers is to exploit the repetitive operation of these devices to update the controller parameters on a cycle-by-cycle basis. The only necessary condition for the application of iterative methods is the existence of measurable variables that allow, albeit indirectly, to evaluate the performance of the system in a given operation. For example, the performance of electromechanical relays has been successfully improved using run-to-run (R2R) algorithms based on measurements of the electrical contacts [9] or the acoustic noise generated when switching [10]. The key idea is to transform the problem into a black-box optimization. For that purpose, the input must be defined by a finite set of parameters and the available measurements must be fed into a cost function that evaluates the performance of the system. The relation between the parameters that define the input signal and the cost is obviously unknown, since it depends both on the system dynamics—which may not be perfectly modeled—and on the possible disturbances and randomness inherent to each operation. For that reason, it is advantageous to use direct search optimization methods to solve the problem, since these do not require information about the gradient and can even be used to optimize discontinuous cost functions. The so-called pattern search algorithms [11] are of particular interest because of their very low computational cost. Bayesian optimization algorithms are also a suitable solution when more computational resources are available [12]. Although iterative algorithms have proven to be a successful solution for improving the performance of electromechan- ical devices, their main disadvantage is that they require a large number of switching operations until convergence is reached. To address this problem, in this paper we present two methods aimed at improving the convergence rate of this class of control algorithms, both based on the sensitivity of the control law with respect to the parameters. Although sensitivity is a well-established and widely used concept in areas such as system identifiability [13] and parameter estimation [14], to the best of our knowledge, it has not yet been applied to improve iterative control algorithms. Specifically, we have applied this methodology to improve the convergence of a recently published adaptive R2R control algorithm for electromechanical relays [15]. The paper is organized as follows. Sections II and III present, respectively, the dynamical model used to design the controller and the R2R algorithm. Section IV explains the sensitivity-based approach to improve convergence, including the two specific methods that we propose. Section V contains simulation results that demonstrate the functionality of our two proposals. Finally, the conclusions are discussed in Section VI. II. SYSTEM DYNAMICS In this section, we provide a brief description of the dynamics of the system to control. For a more detailed explanation, readers are referred to our works [10] and [15]. The electromechanical switching devices under study are single-coil reluctance actuators, consisting of a fixed core wound with a current-carrying coil and a movable iron core. When current flows through the coil, the fixed core becomes magnetized, attracting the movable core. In addition to the magnetic force, the motion of this component is also affected by passive elastic forces, which can generally be modeled as ideal springs. The dynamics of the system can be described using a state-space model, with the position z, velocity v, and magnetic flux linkage λas state variables. The coil voltage, u, is the input. The corresponding state equations are ˙z=v, (1) ˙v=1 m−ks(z−zs)−1 2λ2∂R ∂z ,(2) ˙ λ=−R λ R(z, λ) + u, (3) where mis the moving mass, ksis the spring stiffness, zsis the spring resting position, Ris the coil resistance, and Ris an auxiliary function based on the magnetic reluctance concept. To account for the phenomena of magnetic saturation and flux fringing in the model, Ris defined as R(z, λ) = κ1 1− |λ|/κ2 +κ3+κ4z 1 + κ5zlog(κ6/z),(4) where κ1,κ2,κ3,κ4,κ5, and κ6are positive constants. Overall, the system dynamics depends on q= 9 uncertain parameters, which can be grouped in the parameter vector p. p= [ kszsm κ1κ2κ3κ4κ5κ6]|(5) Note that the resistance Ris treated independently as a parameter without uncertainty, as it can be precisely measured. zd(t)Feedforward controller System Hold Cost Run-to-run adaptation law u[n] ff (t) = uff (t, θ[n]) y[n](t) J[n] θ[n+1] θ[n] Fig. 1. Control diagram. The superscript [n]denotes the variables of the nth operation. The feedforward block computes uff from the parameter vector θand the desired trajectory zd. The adaptation law updates θonce per operation using the cost J, which is derived from the measurable output y. The model (1)–(3) exhibits differential flatness when the position zis regarded as its output. This property allows us to express the state variables and the input as functions of the flat output and its derivatives, as well as the parameter vector p. In particular, the input ucan be expressed as u=fu(z, ˙z, ¨z, ... z , p) = RR(z, λ)λ+˙ λ, (6) where λand ˙ λare given by λ=fλ(z, ¨z, p) = s−2ks(z−zs) + m¨z ∂R ∂z ,(7) ˙ λ=f0 λ(z, ˙z, ¨z, ... z , p) = −ks˙z−m... z−1 2λ2∂2R ∂z2˙z λ∂R ∂z .(8) III. RUN-TO-RUN FEEDFORWARD CONTROL The R2R control strategy presented in [15] and used as the basis for this work is schematized in Fig. 1. It has two main components: a feedforward controller, based on differential flatness theory, and an iterative adaptation law for the controller parameters. A. Feedforward control law The feedforward control law is based on the input function derived using the flatness property, as expressed in (6), which takes a position signal, its derivatives, and the model parameters as inputs. To use this expression as a feedforward control law, a desired position trajectory zdmust firstly be designed. Considering that the objective is to achieve soft landing, zdcan be designed as a 5th-degree polynomial with the following boundary conditions: zd(t0) = z0,˙zd(t0) = 0,¨zd(t0) = 0, zd(tf) = zf,˙zd(tf) = 0,¨zd(tf) = 0,(9) where t0and tfare the desired initial and final times of the switching operation, and z0and zfare the desired initial and final positions, which correspond to the mechanical limits of the motion of the movable core. Once the desired trajectory is set, the feedforward control law can be expressed as a function of a set of parameters with physical interpretation. For reasons related to the adaptation law, each parameter is calculated as the product of its constant nominal value and a dimensionless control parameter. That is, the feedforward control law, uff , is defined as uff (t, θ) = fu(zd(t),˙zd(t),¨zd(t),... zd(t), p∗θ),(10) where θ∈Rqis the vector of control parameters, p∗∈Rq is the nominal parameter vector, and denotes the elementwise—or Hadamard—product. B. Adaptation law One of the drawbacks of feedforward controllers is their dependence on the accuracy of the system model and parameter identification. To mitigate this problem, a run-to-run adaptation law is incorporated. This law acts as an online black-box optimization algorithm that iteratively updates the parameter vector θof the feedforward control law. The optimization objective is to minimize a cost function J calculated from an output variable of the device, namely the absolute value of the impact velocity, denoted as vc. J=|vc|(11) In this context, the chosen optimization method is the pattern search algorithm [11], a derivative-free direct search numerical optimization method. It works by repeatedly applying a specific pattern in a multidimensional search space. Specifically, each evaluated point differs from the central point of the mesh—the best solution so far—in only one of the coordinates. The convergence rate is a key consideration, as each evaluation corresponds to an operation on the real system. To improve this rate, a sensible strategy seems to be to eliminate nearly insensitive coordinates by means of a parametric reduction. This approach is based on the observation that movements in these coordinates result in nearly identical feedforward signals and therefore nearly identical cost values. Thus, avoiding unnecessary exploration in these directions may yield potential benefits. IV. FASTER R2R FEEDFORWARD CONTROL Our proposal in this paper is to improve the convergence rate of the previously described R2R control algorithm by means of a dimensional reduction of the search space. The underlying idea of this approach is that some parameters may have a greater effect than others on the feedforward control law and, as a consequence, on the cost Jto minimize. Thus, by avoiding movements in the directions that have less impact, the algorithm should be able to find close-to-optimal solutions with fewer operations on the real system. In particular, we propose two different reduction methods based on the sensitivity, S(t, θ)∈R1×q, of the feedforward control law (10) to the vector θof control parameters. This sensitivity is given by S(t, θ) = ∂uff (t, θ) ∂θ .(12) Note that, by computing the sensitivity with respect to these adimensional parameters, the elements of S(t, θ)can be compared directly. Otherwise, if the physical parameters had different scales or magnitudes, a normalization would be mandatory. The sensitivity can be used to compute the Fisher information matrix, which is defined as F(θ) = Ztf t0 [S(t, θ)|S(t, θ)] dτ. (13) Note that, by construction, the Fisher matrix F(θ)∈Rq×q is symmetric and positive semidefinite. The first reduction method consists in finding the parameters in θwith the greatest influence on the feedforward control term. Since the sensitivity is a time-dependent signal, let us define the integral-square sensitivity, SIS, as the vector SIS(θ) = Ztf t0 [S(τ, θ)S(τ, θ)] dτ. (14) The components of SIS(θ)∈R1×qquantify the influence of each element of θon uff . Therefore, they can be used to prioritize the parameters according to their relevance and hence to discard the least relevant ones. In this work we assume that the R2R search is conducted in the vicinity of the nominal value of the parameters, i.e., around θ∗= 1q, where 1qis the vector of ones of size q. Thus, the choice of which parameters to optimize and which to keep fixed is determined by SIS(θ∗). The above method provides a way to eliminate the least influential parameters of θ, but it is still possible that two or more parameters classified as highly relevant are somewhat correlated, i.e., have very similar influences on uff . To work around this, the second reduction method is based on finding an alternative orthogonal coordinate system in which the components are ordered according to their influence on the feedforward term. For that, let us define D(θ)as the integralsquare deviation of uff with respect to the nominal input. D(θ) = 1 2Ztf t0uff (τ, θ)−uff (τ, θ∗)2dτ(15) An approximation of D(θ)can be obtained using the secondorder Taylor expansion around the nominal parameter vector, D(θ)≈D(θ∗) + Dθ(θ∗)δθ +1 2δθ|Dθθ(θ∗)δθ, (16) where δθ =θ−θ∗, Dθ(θ) = ∂D ∂θ , Dθθ(θ) = ∂2D ∂θ2. It is not difficult to see that, while both D(θ∗)and Dθ(θ∗) are equal to zero, the Hessian matrix of Dat θ∗is equal to the Fisher matrix evaluated at the same point, i.e., Dθθ(θ∗) = F(θ∗). Therefore, the integral-square deviation of uff with respect to the nominal input is approximately given by the quadratic form D(θ)≈1 2δθ|F(θ∗)δθ. (17) Now, the goal is to find a linear change of variables to an orthogonal coordinate system in which each component affects Dindependently. Since F(θ∗)is symmetric and positive semidefinite, the eigendecomposition and the singular value decomposition coincide. That is, the Fisher matrix can be expressed as F(θ∗) = VΛV|= q X i=1 viλiv| i,(18) where V∈Rq×qis an orthogonal matrix whose columns, vi∈Rq, are the eigenvectors (or singular vectors) of F(θ∗) and Λ∈Rq×qis a diagonal matrix with the corresponding eigenvalues (or singular values), λi∈R, in the diagonal elements. Assuming that the eigenvalues are sorted by value, λ1≥λ2≥ · · · ≥ λq≥0,(19) the Fisher matrix could be approximated by F(θ∗)≈˜ F(θ∗) = ˜ V˜ Λ˜ V|= r X i=1 viλiv| i,(20) where ris the order of the approximation (r≤q) and ˜ V∈Rq×rand ˜ Λ∈Rr×rcontain the first reigenvectors and eigenvalues, respectively. Using this reduced-order approximation, Dcan be further approximated as D(θ)≈1 2δθ|˜ F(θ∗)δθ =1 2(ϕ−ϕ∗)|˜ Λ (ϕ−ϕ∗),(21) where ϕ∈Rris the alternative orthogonal reduced-order parameter vector, given by ϕ=˜ V|δθ +ϕ∗=˜ V|(θ−θ∗) + ϕ∗,(22) and ϕ∗is its nominal value, which can be chosen arbitrarily. By selecting ϕ∗=˜ V|θ∗, the change of variables simplifies into ϕ=˜ V|θ⇐⇒ θ=˜ V ϕ. (23) The R2R search in this second method is then performed in the reduced space defined by the alternative parameterization ϕ. As in the first method, the search is assumed to be performed in the vicinity of the nominal parameter vector, so the change of basis, which is fixed, is based on the decomposition of F(θ∗). V. SIMULATION RESULTS In this section, we implement the proposed parameter reduction techniques to obtain different parameterizations of the feedforward term. We evaluate the control strategy for each parameterization through simulated experiments. For all the experiments, the feedforward signal uff is initialized based on the nominal parameter values specified in Table I. For the sake of brevity, we focus on the closing operation, i.e., z0and zfare set as the upper and lower position limits of the stroke of the movable core. The control process of the opening operation is completely equivalent. A. Dimensional reduction of the search space As first step, we have computed the integral-square sensitivities, SIS, of the feedforward term to both the original control parameters and those resulting from orthogonalization. These sensitivities, which are displayed in Fig. 2, provide insights into potential parametric reduction cases. As previously presented, the first reduction technique requires computing the sensitivity with respect to each original parameter θi(the subscript idenotes the ith vector component). These sensitivities, which are visualized in Fig. 2a, indicate two main possible reduction cases: eliminate from the search only θ4and θ6, which are by far the least TABLE I NOMINAL PARAMETER VALUES ks55 N/mκ51320 m−1 zs0.015 m κ69.73 ·10−3m m1.6·10−3kg R50 Ω κ11.35 H−1z010−3m κ20.0229 Wb zf0 κ33.88 H−1t00 κ47.67 ·104H−1/mtf3.5·10−3s (a) Original parameterization. (b) Orthogonal parameterization. Fig. 2. Integral-square sensitivities of the feedforward control law with respect to the control parameters. influential parameters; or eliminate θ7,θ8, and θ9as well, because their SIS values are also relatively small. The second reduction technique involves a change of basis resulting in a new parameterization based on the orthogonal parameters. The sensitivities to these new parameters have been computed and represented in Fig. 2b. Note that these sensitivity values, which are the eigenvalues of the Fisher matrix F(θ∗), have been represented in order from the most to the least influential in the feedforward control law. Based on the visualized values, a new possible reduction case is suggested: eliminate from the search all but ϕ1and ϕ2, as they are the most influential ones by a wide margin. In summary, we have identified three promising reducedorder parameter vectors, with dimensions r= 7 and 4for the original parameters, and r= 2 for the orthogonal parameters. Nonetheless, in order to ensure a comprehensive comparison, in the following we evaluate equivalent parametric reductions, with r= 7,4, and 2, to both sets of parameters. We also consider a baseline scenario in which no reduction is performed, i.e., r=q= 9. All studied reduction cases are summarized in Table II. B. Control results The controller performance with the proposed reductions has been evaluated through simulations using the dynamical model presented in Section II. We assume that the controlled dynamic system is governed by these equations. However, considering that the parameters are never perfectly known in practice, it is assumed that there is some uncertainty in their values. In particular, each component of the model parameter vector pis randomly and independently perturbed between 95 % and 105 % of its nominal value (see Table I). Note TABLE II REDUCTION CASES Case rFree parameters Fixed parameters A9θ1, . . . , θ9 B7θ1, θ2, θ3, θ5, θ7, θ8, θ9θ4, θ6 C4θ1, θ2, θ3, θ5θ4, θ6, θ7, θ8, θ9 D2θ2, θ3θ1, θ4, θ5, θ6, θ7, θ8, θ9 E7ϕ1, . . . , ϕ7ϕ8, ϕ9 F4ϕ1, . . . , ϕ4ϕ5, . . . , ϕ9 G2ϕ1, ϕ2ϕ3,..., ϕ9 that the electrical resistance Ris not perturbed because it is typically measured with high precision. For each reduction case, 10 000 different trials have been simulated, and in each trial the control algorithm is run for 300 switching operations. The feedforward control law, on the other hand, starts all trials using the nominal value, i.e., θ[1] =θ∗= 1q. The control results for each case are summarized in Fig. 3, where each graph represents the obtained distribution of costs, J=|vc|, with respect to the switching operation, n. To show the effectiveness of each control scenario, the graphs also display the cost of a conventional switching operation, namely with a 30 V constant activation. For the first simulated case, serving as a baseline, no parameter reduction is applied (case A from Table II). Note that this scenario corresponds to the direct application of the control algorithm presented in [15]. As shown in the control results (see Fig. 3a), the disparity between the system parameters (randomly perturbed) and the initial feedforward parameters results in a large variability in the impact velocities of the first switching operation (n= 1). Then, due to the adaptation law, the control performance improves as the number of iterations increases. However, it is evident that the control convergence is quite slow, as it requires more than 300 operations to reach stable values of the cost. The main reason for this is the large number of parameters modified by the adaptation law. Subsequently, the parameter reductions are applied with the aim of improving the controller convergence. In the next simulated cases, the reduced original parameter vectors (with dimensions r= 7,4, and 2) are used for the adaptation law. The control results are displayed in Figs. 3b, 3c, and 3d, respectively. As can be seen, the number of iterations required to improve the performance is reduced. In particular, with no parameter reduction, the control needs to perform 203 operations to halve the cost of the uncontrolled scenario in 90 % of the trials. In contrast, with this type of parameter reduction, the control strategy only requires 153, 83, and 42 operations in each subsequent case. Furthermore, the results with the reduced orthogonal parameters are visualized in Figs. 3e, 3f, and 3g for r= 7, 4, and 2, respectively. These also demonstrate a significant improvement in the control convergence rate with respect to the baseline. However, when each case is compared with its equivalent case (i.e., with the same r) using the (a) Case A: No reduction (r= 9). (b) Case B: Original parameters (r= 7). (c) Case C: Original parameters (r= 4). (d) Case D: Original parameters (r= 2). (e) Case E: Orthogonal parameters (r= 7). (f) Case F: Orthogonal parameters (r= 4). (g) Case G: Orthogonal parameters (r= 2). Fig. 3. Cost values with respect to the number of switching operations. Each graph shows the median (P50) and the 10th and 90th percentiles (P10 and P90, respectively) of the distribution of values obtained for the 10 000 simulated experiments. The cost without control is also represented. original parameterization, a significant similarity is observed, making it challenging to discern the superior performer. For instance, concerning r= 7, case E (with the orthogonal parameterization) presents better results on average than case B (using the original parameters). However, it also exhibits a higher variability. This suggests that the poor performance of case E may be due to the negligible effect of the last included orthogonal parameters, ϕ5to ϕ7, on the feedforward control law. In contrast, for r= 4 and r= 7 the results with both reduction techniques are very similar. To facilitate comparisons across all the studied cases, we have computed the integrated (i.e., cumulative) cost of each trial, denoted as I, which provides a measure of the overall performance of each control strategy for any iteration count. I[n]= n X i=1 J[n](24) The average values of these integrated costs across all 10 000 trials are represented in Fig. 4. For any count of the tested operations, from 1 to 300, it is evident that reducing the dimension rof the search space improves the total cost. This representation also shows that reductions with orthogonal parameters have superior performance for r= 7 and r= 4, although it is slightly less effective for r= 2. Moreover, by examining the Ivalues and their slopes for the last operation (n= 300), the main drawback of extreme parametric reduction (r= 2) becomes apparent: if the number of operations were to increase further, the total control cost would grow faster than in the other cases, because these have converged to worse stable values. These results also demonstrate that, for a larger number of iterations, the best performance is achieved with the orthogonal parameterization of dimension r= 4. In any case, the results show that the control convergence is improved in all the considered scenarios, which emphasizes the potential of the proposed sensitivitybased reduction methodology. VI. CONCLUSIONS In this work, we have analyzed the possibility of improving the convergence rate of an adaptive feedforward R2R control algorithm by means of a sensitivity-based dimensional reduction of the search space. The two strategies presented show great potential in simulation, but the results obtained do not yet allow us to determine whether one of the two is more efficient than the other. Possibly, this will depend on the system to be controlled and the physical parameterization originally chosen for the dynamical model. As future work, we would like to evaluate the possibility of periodically updating the reduced parametric basis. That is, the two proposals of this work use a fixed transformation based on the nominal value of the parameters, but it would be interesting to evaluate whether it is advantageous to redo the calculations every iteration or every few iterations. In addition, we also intend to perform real laboratory tests and verify that the experimental results agree with those observed in simulation. Fig. 4. Integrated costs with respect to the number of iterations. Each line represents the mean values for the 10 000 simulated experiments. Comparison between different reduction methods. REFERENCES [1] J. van Dam, D. Kropl, B. Gysen, and E. Lomonova, “Design of a permanent magnet-biased reluctance valve actuator with integrated eddy current damping,” in Int. Symp. Linear Drives for Ind. Applic. IEEE, 2019, pp. 1–6. [2] A. Di Gaeta, C. I. Hoyos Velasco, and U. Montanaro, “Cycle-by-cycle adaptive force compensation for the soft-landing control of an electromechanical engine valve actuator,” Asian J. Control, vol. 17, no. 5, pp. 1707–1724, Jul. 2015. [3] B. Gergiˇ c and D. Hercog, “Design and implementation of a measurement system for high-speed testing of electromechanical relays,” Measurement, vol. 135, pp. 112–121, Mar. 2019. [4] F. Deschaux, F. Gouaisbaut, and Y. Ariba, “Magnetic force modelling and nonlinear switched control of an electromagnetic actuator,” in IEEE Conf. Decision and Control. IEEE, 2019, pp. 1416–1421. [5] J. Fang, X. Wang, J. Wu, S. Yang, L. Li, X. Gao, and Y. Tian, “Modeling and control of a high speed on/off valve actuator,” Int. J. Automot. Technol., vol. 20, pp. 1221–1236, 2019. [6] M. Benosman and G. M. Atınc¸, “Extremum seeking-based adaptive control for electromagnetic actuators,” Int. J. Control, vol. 88, no. 3, pp. 517–530, 2015. [7] W. Hoffmann, K. Peterson, and A. G. Stefanopoulou, “Iterative learning control for soft landing of electromechanical valve actuator in camless engines,” IEEE Trans. Control Syst. Technol., vol. 11, no. 2, pp. 174–184, Mar. 2003. [8] T. Braun, J. Reuter, and J. Rudolph, “Flatness-based feed-forward control design for solenoid actuators considering eddy currents,” IFACPapersOnLine, vol. 52, no. 15, pp. 567–572, 2019, 8th IFAC Symp. Mechatronic Syst. [9] E. Ramirez-Laboreo, C. Sagues, and S. Llorente, “A new run-to-run approach for reducing contact bounce in electromagnetic switches,” IEEE Trans. Ind. Electron., vol. 64, no. 1, pp. 535–543, Jan. 2017. [10] E. Serrano-Seco, E. Ramirez-Laboreo, E. Moya-Lasheras, and C. Sagues, “An audio-based iterative controller for soft landing of electromechanical relays,” IEEE Trans. Ind. Electron., vol. 70, no. 12, pp. 12 730–12 738, Dec. 2023. [11] R. M. Lewis and V. Torczon, “Pattern search methods for linearly constrained minimization,” SIAM J. Optimization, vol. 10, no. 3, pp. 917–941, 2000. [12] E. Moya-Lasheras, E. Ramirez-Laboreo, and C. Sagues, “A novel algorithm based on Bayesian optimization for run-to-run control of short-stroke reluctance actuators,” in European Control Conf. IEEE, Jun. 2019, pp. 1103–1109. [13] H. Miao, X. Xia, A. S. Perelson, and H. Wu, “On identifiability of nonlinear ode models and applications in viral dynamics,” SIAM review, vol. 53, no. 1, pp. 3–39, 2011. [14] E. Walter, L. Pronzato, and J. Norton, Identification of parametric models from experimental data. Springer, 1997, vol. 1, no. 2. [15] E. Moya-Lasheras, E. Ramirez-Laboreo, and E. Serrano-Seco, “Runto-Run Adaptive Nonlinear Feedforward Control of Electromechanical Switching Devices,” IFAC-PapersOnLine, vol. 56, no. 2, pp. 5358– 5363, 2023, 22nd IFAC World Congr.