Design Methodology for a Magnetic Levitation System Based on a New Multi-Objective Optimization Algorithm

Abstract

Multi-objective (MO) optimization is a developing technique for increasing closed-loop performance and robustness. However, its applications to control engineering mostly concern first or second order approximation models. This article proposes a novel MO algorithm, suitable for the design and control of mechanical systems, which does not require any order reduction techniques. The controller parameters are determined directly from a special type of rapid analysis of simulated transient responses. The case study presented in this article consists of a magnetic levitation system. Certain difficulties such as the nonlinearity identification of the magnetic force and duo magnetic field sensor scheme were addressed. To point out the advantages of using the developed approach, the simulations as well as the experiments performed with the help of the created algorithm were compared to those made with common MO algorithms.

1. Introduction

Any engineering system goes through a design stage. A general rule is to try to resolve as many obstacles as possible during the design stage—this includes controller tuning. Certain circumstances may exist that interfere with online tuning. For example, the plant may have to be put offline. This may cause unnecessary stalls for a production line, which this particular system is a part of, which, depending on the process at hand, could be costly. Another case may be a long loop time—some processes in industry take a matter of days or longer. During the tuning process, it is naturally required to have several step-responses for a better quality of the tune. Finally, it is beneficial to have an idea of a good tune before the final implementation, which may lead to a better optimization overall.

The problem of applying and optimizing a proper controller is one of the central tasks of engineering. Whether this is a classic PID controller or a nonlinear one, the goal is to help the physical process at hand run within the needed boundaries and properties. A revisit of the classic control theory and a re-evaluation of its possibilities with modern computing powers provides for a new controller tuning method and a powerful tool for design optimization. Our method of controller tuning deals with a plant in the design stage. An algorithm was developed which performs a search in the parameter space with all the main transient response characteristics being rapidly estimated by it. It is based on fast computations of the inverse Laplace transform and curve analysis of the continuously generated stochastic Laplace images.

The idea behind this method appeared while working with a magnetic levitation system. A demand for a more efficient tune of a PID controller during the design stage led to the development of the algorithm described in this work. PID controllers are widely used in various industrial applications. The effectiveness of such controllers depends on tuning, which is essentially an optimization problem. To address this task, a number of tuning methods have been developed ever since the appearance of these controllers [1,2]. Since this system is open-loop unstable, experimental methods are difficult to implement, and it can be frustrating to tune it manually without a model. On the one hand, classic control has powerful tools for system design such as root locus, frequency response method, etc. [3,4]. However, the developed algorithm is based on a different approach. In a situation where the developed system’s model is known—fast computations of the inverse Laplace transform are used to provide lots of information about the studied system’s behavior. Whether it is the classical partial fraction expansion [5] or a numerical method [6,7,8] (as long as it is fast enough), combined with the developed algorithm, they make for a unified framework for designing and tuning closed-loop controlled mechanical systems.

The goal of this work was to explore this possibility and present an overview of such an approach. The developed algorithm is suitable for multi-objective optimization. Since a transient response of a mechanical system has a few characteristics that may be important for the designer, controller tuning can be formulated as a multi-objective problem. The developed algorithm is compared with a common multi-objective algorithm NSGA-II using Pareto-optimal front. At the same time, we discuss the challenges of the magnetic levitation system itself such as the identification of the magnetic nonlinearity. A detailed physical background is given in order to obtain the system model. At the end, we show how this approach is a foundation of a larger problem of system design, as this method can be applied to solve an inverse problem of finding the optimum plant parameters for the desired transient response characteristics.

2. Materials and Methods: The Magnetic Levitation System

The magnetic levitation set-up involves a well-known phenomenon often used for control system studies. Recently, it has become important in a very wide range of industrial applications where magnetic suspension techniques can be profitably applied. The best known ones are high-speed ground transportation [9,10] and high-speed bearings with reduced noise and friction [11,12].

Levitation can in general be achieved in two ways. The first one is using AC in a primary coil. As a result a current is induced in a secondary coil, which is repelled from the primary one [13]. The height of the levitation can be controlled by amplitude and/or frequency of the AC. The second method is using DC current in the primary coil and a permanent magnet or a piece of ferromagnetic material as the levitated object. If a permanent magnet is being used, both attraction and repulsion are feasible. If however a ferromagnetic material is being used for levitation, only attraction is possible. From a control theory point of view it makes a difference whether attraction or repulsion is used in a magnetic levitation project (see Figure 1).

In the case of repulsion, the open loop system is stable. If, however, attraction is being used, the system is open loop unstable and is useless without a proper closed-loop design. In order to study magnetic levitation control, various experimental setups have been used  [15,16] but the most commonly used experimental setup is the one shown in Figure 2 [17,18].

The main goal of such systems is to make the permanent magnet levitate at a desired height. Various approaches known from control theory have been used for this purpose (root locus [19,20], state space [21], disturbance rejection control [22], fuzzy control [23], sliding mode control [24], fuzzy sliding mode control [25], robust control [26], neural network control [27], various nonlinear approaches [28,29]). The authors of [30] deal with dynamical uncertainties and exterior perturbations in a magnetic levitation system using a real-time prescribed performance control. This allows for chattering reduction and faster convergence to the equilibrium point. In [31], an analytical method using Lagrange equations for the analysis of magnetic levitation (MagLev) systems is proposed. This provides for an interesting MagLev model which distinguishes the primary and induced currents and also the equilibrium height of the levitating object on the input voltage through the mutual inductance of the system.

Before we discuss the designed controller tuning and system design method, let us describe the constructed magnetic levitation system. Let us start with the development of the system’s model. To do that, we use the classic block diagram representation shown in Figure 3.

The $z_d$ is the desired position of the permanent magnet (input signal), while z is the actual position (output signal). The $G_c,G_a,G_o$ and $G_s$ are the transfer functions of the corresponding parts of the experimental setup. Application of block diagrams leads to very descriptive relationships, especially in the case of automatic feedback systems. To acquire the model of the system we need to determine all of the expressions behind these blocks.

2.1. Object

By far the most difficult element to model is the object. Its input signal is the voltage $U_c$ on the coil and the output signal is the position of the magnet z. First, we should analyze the force acting upon a permanent magnet in a magnetic field.

2.1.1. Determination of the Magnetic Force

To determine the magnetic force $F_m$, we used the experimental setup shown in Figure 4. An electromagnetic coil is a length of wire wound in a joined sequence of concentric rings through which an electric current i flows. The magnetic field $\mathbf{B}$ (or its component $B_z$) of a single ring can be obtained by the application of the Biot–Savart law,

(1)$B_z=\frac{{\mu }_0}{2\pi }i\frac{\pi R^2}{{(z^2+R^2)}^{3/2}},{\frac{\partial B_z}{\partial z}}_{coil}=-\frac{3}{2}·\frac{{\mu }_{0}i{R}^{2}z}{{(z^2+R^2)}^{5/2}},$

where ${\mu }_0=4\pi ·{10}^{-7}$ H/m is the magnetic permeability of vacuum and R is the radius of the coil’s winding.

A permanent magnet can be modeled as a collection of many microscopic current loops (magnetic dipoles). The net effect of these small current loops is a surface current $i_m,$ which is called the Amperian current [32]. Let the current loop (magnetic dipole) have a magnetic moment of $\mathbf{\mu }=\left({\mu }_x,{\mu }_y,{\mu }_z\right)$ and be in a uniform magnetic field $\mathbf{B}=\left(B_x,B_y,B_z\right).$ If the loop is small enough, then the torque acting on such a loop is given by a simple expression [33]:

(2)$\mathit{\tau }=\mathit{\mu }\times \mathbf{B};$

therefore, the force $\mathbf{F}$ acting on a magnetic dipole in this field is the gradient of the potential energy associated with this torque:

(3)$\mathbf{F}=-\nabla U_B=-\nabla \left(-\mathit{\mu }·\mathbf{B}\right)=\nabla \left(\mathit{\mu }·\mathbf{B}\right).$

As magnetic dipole moment $\mathbf{\mu }$ has only a vertical component (${\mu }_z$ in this case), which is constant, we can write the vertical component of the force $F_z$ as follows:

(4)$F_z={\mu }_z\frac{\partial B_z}{\partial z}.$

The magnet was attached to the plastic pedestal using adhesive tape. The pedestal was made with a screw so that it is easy to put the magnet at a desired height z. The pedestal was put on scales. When there is an attractive force between the electromagnet and the permanent magnet, the reading on the scales is lower. In this way, a matrix of possible heights z (measured from the lower end of the electromagnet) and electrical currents i through the coil was formed. The results are given in Table 1.

The relationship between these variables can be obtained using any available curve fitting tool. The data from Table 1 must be imported and then fitted using a custom function f in the form of $F_m=f(i,z)$, where i and z are independent variables of current and coordinate and $F_m$ is the dependent variable of magnetic force. In many papers dealing with the magnetic levitation, the force of the coil is approximated by the formula $F_m=const·i/z^2$. This yields acceptable results in some cases but in this research we paid extra attention to the accuracy of the acquired expression.

For a real magnetic coil, it is often quite difficult to acquire an explicit expression of the magnetic field. One could be surprised by the scale of difficulties appearing on this path. For objects with an infinite dimension, explicit formulas usually do exist since it is possible to perform a limit passage. Another possibility is the geometry of the given current contour having a symmetry axis (like the solenoid when the field is calculated at some point on the coaxial line of the magnet). One of the obstacles to consider is the fact that the majority of the formulas are for a magnetic force applied to a material point. In the case of a magnetic force applied to a permanent magnet things are even more complicated due to the shapes of the interacting magnetic objects as effects such as mutual inductance have to be taken into account. In the case of non-simplest volumetric bodies it is often impossible to integrate the expression (1) analytically. Numerous approaches exist to simplify this process such as Maxwell’s Method [34]. This is why, in engineering, numerical approaches such as the finite elements method or the boundary integral equations method are often implemented. However, in our research we require an explicit expression for the force applied to the permanent magnet.

It is safe to assume that a model with an explicit formula representing a closer geometry to the original set up would provide a better fit for the given measured data; we show that in Table 2. Here, we establish that one of the best fits for our data is, in fact, the expression using (4):

(5)$F_z=-a·\frac{3{\mu }_z}{2}·\frac{{\mu }_{0}i{R}^{2}z}{{(z^2+R^2)}^{5/2}},$

where ${\mu }_0=4\pi ·{10}^{-7}$ H / m is the magnetic permeability of vacuum, ${\mu }_z$ is a magnetic dipole moment of the permanent magnet and a is a constant related to the length of the coil L and the winding turns per unit length n. The length of the coil is $L=54\mathrm{mm}$, while the mean radius of the coil’s winding is $R=37\mathrm{mm}$. The permanent magnet has a form of a cylinder with 4 mm radius and 5 mm height. For a neodymium magnet of this size, the vertical component of the magnetic moment can be estimated to have a value of $|{\mu }_z|=0.49\mathrm{A}·{\mathrm{m}}^2$. As a means of comparison we chose the commonly used functions:

(6)$F_z=a\frac{i}{z},F_z=a\frac{i}{z^2},F_z=a\frac{i}{z^3}.$

Using curve fitting software, we determined the value of parameter a and also the plausibility of the fit as a sum of squared errors (SSE). The results are summarized in Table 2.

The fit using expression (5) has a sum of squared errors of $5.6\times {10}^{-5}$, which is by two orders of magnitude lower than the most used fit $ax/y^2$. Therefore, expression (5) provides for a much better approximation.

2.1.2. Experimental Determination of Coil Resistance and Inductance

The electromagnet used in our experiment was a reel of PU enamelled, unjacketed copper wire. Its length was $230\mathrm{m}$ and the cross sectional area was $0.246{\mathrm{mm}}^2$. The resistance of the coil $R_c$ can easily be measured using a multimeter. In the case of our coil it was $R_c=16.3\Omega$. In order to determine the inductance of the coil $L_c$, a resistor $R_r=468\Omega$ was added in series to the coil as shown schematically in Figure 5.

The ratio between output voltage $U_{out}\left(t\right)$ and input voltage $U_{in}\left(t\right)$ can be obtained by a simple voltage divider rule,

(7)$\frac{U_{out}\left(j\omega \right)}{U_{in}\left(j\omega \right)}=\frac{Z_c}{Z_r+Z_c}=\frac{R_c+j\omega L_c}{R_r+R_c+j\omega L_c}.$

The ratio of amplitudes of the output and input voltages is, therefore, given by:

(8)$\frac{U_{out}}{U_{in}}=\frac{|U_{out}\left(j\omega \right)|}{|U_{in}\left(j\omega \right)|}=\frac{\sqrt{R_c^2+{\left(\omega L_c\right)}^2}}{\sqrt{{(R_r+R_c)}^2+{\left(\omega L_c\right)}^2}}.$

By a simple algebraic manipulation, we get the following result:

(9)$L_c=\frac{1}{\omega }\sqrt{\frac{U_{out}^2{(R_c+R_r)}^2-U_{in}^2{R}_c^2}{U_{in}^2-U_{out}^2}}.$

The actual waveforms for both $U_{in}\left(t\right)$ and $U_{out}\left(t\right)$ are shown in Figure 6.

For the signal $U_{in}=5.12\mathrm{V}$ and $\omega =6280\mathrm{rad}/\mathrm{s}$ we get the $U_{out}=2.87\mathrm{V}$. Finally, using Equation (9) we get $L_c=52.1\mathrm{mH}$.

2.1.3. The Transfer Function of the Object

With all the parameter values known, we can now create the model of the object, the input of which is the voltage while the output is the current $i\left(t\right)$.

(10)$G_{coil}\left(s\right)=\frac{\widehat{I}\left(s\right)}{\widehat{U_c}\left(s\right)}=\frac{1}{Ls+R_c}.$

If we neglect the viscous resistive forces of air, the only forces acting on the magnet are the gravity and the magnetic force. The movement in the vertical direction (z-axis) can then be described by the following differential equation:

(11)$m\frac{d^{2}z}{d{t}^2}=F_z-mg,$

where $m=3\mathrm{grams}$ is the mass of the magnet, g is the gravitational acceleration and the force $F_z$ is given by the expression (5). Let $\widehat{z}\left(t\right)$ be the relative change of coordinate z from the initial state $z_0$:

(12)$z\left(t\right)=z_0+\widehat{z}\left(t\right).$

Zero initial value $\widehat{z}(+0)=0$: satisfies the rule for differentiation of originals for the Laplace Transform [5]. It is necessary to linearize the force $F_m$ in Equation (11). We select the operating point to be the equilibrium state at $z_0=-25\mathrm{mm}$. Using (5) to calculate the initial current from equation $F_z=mg$ we get $i_0=0.23\mathrm{A}$. The initial voltage is $U_{c0}=i_0{R}_c=3.67\mathrm{V}$.

We expand the function $F_z$ in a Taylor series at the point $(i_0,z_0)$, where $\widehat{i}=0,\widehat{z}=0$. Since the operating point is chosen in a way that $mg=F_z(i_0,z_0)$ we rewrite the Equation (11) as

(13)$m\frac{d^2\widehat{z}}{d{t}^2}=\frac{\partial F_z}{\partial z}(i_0,z_0)·\widehat{z}+\frac{\partial F_z}{\partial i}(i_0,z_0)·\widehat{i},$

where

(14)$\frac{\partial F_z}{\partial z}(i_0,z_0)=-\frac{3{\mu }_z{R}^2{i}_0}{2}\frac{R^2-4z_0^2}{{(z_0^2+R^2)}^{7/2}}$

and

(15)$\frac{\partial F_z}{\partial i}(i_0,z_0)=-\frac{3{\mu }_z{R}^2}{2}\frac{z_0}{{(z_0^2+R^2)}^{5/2}}.$

The Laplace Transform of the Equation (13) is

$s^2\widehat{Z}\left(s\right)=a\widehat{Z}\left(s\right)+b\widehat{I}\left(s\right),$

where

$a=\frac{1}{m}\frac{\partial F_z}{\partial z}(i_0,z_0)=-783.3,b=\frac{1}{m}\frac{\partial F_z}{\partial i}(i_0,z_0)=39.9.$

Note that, with this system and the direction of the Z-axis, a and b must be greater than $0.$ The transfer function of Equation (13) is, therefore:

$G_{magnet}\left(s\right)=\frac{\widehat{Z}\left(s\right)}{\widehat{I}\left(s\right)}=\frac{b}{s^2-a}.$

Combined with the expression (10) the resulting transfer function of our magnetic levitation system is

(16)$G_o\left(s\right)=\frac{b/L_c}{s^3+(R_c/L_c)s^2-as-a{R}_c/L_c}=\frac{766}{s^3+313s^2-783s-2.45\times {10}^5}.$

2.2. Actuator

As shown in Figure 3, an actuator is an element between the controller and the object. In our case, its primary role was to supply proper voltage and current for the electromagnet. It is shown schematically in Figure 7.

It was made of an optocoupler (LED—phototransistor pair and a MOSFET driver) and a MOSFET. The MOSFET was in series connection with the electromagnet. A snubber diode was added in order to suppress the transients, which appeared due to the pulse width modulated input signal $U_{PWM}$ (which is the output voltage of the controller shown in Figure 7). The maximum value of $U_{PWM}$ is equal to $3.33\mathrm{V}$. The voltage $U_{cc}$ was selected to be $16.2\mathrm{V}$, so that the maximum current through the electromagnet is $1\mathrm{A}$. The transfer function of the actuator is, therefore, equal to:

(17)$G_a=\frac{16.2\mathrm{V}}{3.33\mathrm{V}}=4.86.$

2.3. Sensor

The role of the sensor is to detect the actual position of the permanent magnet. In general it can be detected in two ways as shown in Figure 8.

The first and the most common one is the application of a photo-emitter and a photo-receiver [17,25,28,35,36,37] or some similar arrangement based on optical means [22,38] The second one is the application of one or two Hall sensors [39,40]. Sometimes an inductive sensor is used [41,42]. In our case, a pair of SS49E Hall sensors have been used. The idea behind using a Hall sensor is to detect the magnetic field of the permanent magnet. When it moves closer to the electromagnet, an increase in the magnetic field can be detected. The problem is, however, that the electromagnet generates a magnetic field of its own as well. When the current through the electromagnet is changed due to the varying control signal, the sensor cannot tell if the magnetic field changed due to the changed current through the electromagnet or due to the movement of the permanent magnet. If a pair of such sensors is used (one at the upper and one at the lower end of the electromagnet), their output signals can be subtracted and the resulting signal is due to the magnetic field of the permanent magnet alone. In order to get this signal the circuit shown in Figure 9 is used.

It is composed of two operational amplifiers. The first one is used for subtraction and the second one for amplification. The circuit output for various distances of the permanent magnet is given in Table 3.

The dependence between the voltage $U_s$ and the permanent magnet position can then be obtained in a similar way as in the case of the magnetic force $F_m$. Using a curve fitting software, we get a reasonably good approximation with:

(18)$U_s=f_s\left(z\right)=\frac{-1.707\times {10}^5}{z^4}.$

At this stage, we could also make linearization and get the $G_s$ denoted in Figure 3, but it is easier to include $z=f^{-1}\left(U_s\right)$ in the controller.

2.4. Controller

The role of the controller is played by the Arduiuno Due microcontroller. In reality, the block diagram of the whole system shown in Figure 3 should be modified as shown in Figure 10.

The microcontroller has two inputs. One is the sensor voltage $U_s$, which can be used to obtain the permanent magnet position z using

(19)$z=f_s^{-1}\left(U_s\right)=-\sqrt[4]{\frac{1.707·{10}^5}{U_s}}.$

The other one is the desired position of the permanent magnet $z_d$, which can be entered into a microcontroller via a Serial Monitor (serial communication with a PC). With the subtraction of z from $z_d$ we get the error signal. The error signal is the input for $G_c$ (block Contr. Alg. in Figure 3 and Figure 10). The control algorithm is the common PID algorithm. The control signal (output of the controller) is a pulse width modulated signal.

3. Materials and Methods: The Developed Algorithm

3.1. Problem Statement

Providing the underlying calculation method is fast enough, one can look at the tuning problem from a different angle. While the classical methods do exist and provide the user with satisfying results, sometimes they seem a little bit locked on to their sequence of actions. It is always a good idea to let the simulation “run free” in terms of possible perturbations, irregularities and parameter values. Real systems tend to yield slightly different results to simulations. This raises an important question—what is the “optimal” controller tuning? Should we stop when have reached acceptable transient response characteristic or is it beneficial to keep going to explore the system’s behavior over an area of controller parameters? By area we mean certain intervals within the parameter space (rectangles on a plane, certain cubic areas in 3-dimensions parameter space).

In the literature, the transient response analysis mostly comes down to analyzing the known characteristics of a first- or second-order-like mechanical system [43,44,45,46,47]. Usually, systems of higher order are approximated to it by known techniques. The convenience of the second-order-like system approximation, the simplicity of the action–result methodology of adjusting loop gains led to other techniques of step-response analysis being overlooked. Higher order systems do not have such an intuitive correlation between transient response parameters and controller gains. Therefore, a second-order approximation is commonly used.

In this section, we describe a new optimization algorithm and test it in a case of the magnetic levitation of a small cylinder with a PID controller as means of keeping it at a desired height. PID controllers are still one of the most commonly used controllers in industrial applications [1]. The idea behind them is intuitive—with the help of examples, one can understand the core principles without referring to the Laplace Transform. There are numerous ways of tuning a PID controller. Many of these methods involve the step-response method either by using a process model or an experiment. The output is measured or calculated as a function of time. By analyzing it, a new set of controller parameters is chosen. Many strategies exist in properly applying a PID controller. Ref. [48] suggests a heuristic algorithm using wavelets for online tuning of a gain adapted PID-controlled linear actuator. Permanent magnets are implemented as excitation with the aim of soft landing which increases reliable functionality and component life. In [49], another technique is utilized to achieve tracking and soft landing for electromagnetic actuators. Pre-action is employed to enable the system to avoid power saturation. For constrained processes, it was demonstrated that a PID with anti-windup is able to provide similar or even better results than model predictive control when certain solutions are considered [43]. Conditions on nonlinearity and uncertainty are addressed in [50] so that a high order affine-nonlinear system under an extended PID controller can be semi-globally stabilized with a fast rate of regulation error convergence.

First, let us unite the actuator and object into a single plant block. The resulting transfer function of the plant would be $G_{ao}\left(s\right)=G_a\left(s\right)·G_o\left(s\right),$ and

(20)$G_{ao}\left(s\right)=\frac{3723}{s^3+313s^2-783s-2.45\times {10}^5}.$

Due to negative coefficients of the polynomial in the denominator the system is unstable, hence a controller is needed. As seen in Figure 10, since $G_s\left(s\right)·G_s^{-1}\left(s\right)=1$, we have a unity feedback control system. The transfer function of the PID controller in a parallel form is $G_c\left(s\right)=K_p+K_i/s+K_d·s,$ where $K_p,K_i$ and $K_d$ are positive real values. The resulting transfer function of the controlled magnetic levitation system $T\left(s\right)$ ends up being

(21)$\frac{A(K_d{s}^2+K_{p}s+K_i)}{a_3{s}^4+a_2{s}^3+(a_1+A{K}_d)s^2+(a_0+A{K}_p)s+A{K}_i},$

where $A=3723,a_3=1,a_2=312.9,a_1=-783.3,a_0=-2.45·{10}^5$. This transfer function is going to be the subject for testing our algorithm.

In order to obtain a step response of our system, first we multiply the transfer function (21) by $1/s$ and then perform the inverse transformation. Since the function (21) is in the form of a polynomial divided by a polynomial of degree lesser than the one in the nominator, then the inverse Laplace transformation of such a function is given by a partial fraction expansion [5]. Therefore, the explicit expression of the signal $f\left(t\right)$ and its derivative $f^{\prime }\left(t\right)$ at any given point in time is presented.

We are going to show how, for this magnetic levitation system represented as a transfer function, our algorithm will provide numerous stable solutions while recording all the necessary step-response characteristics. An array with such data is created in the process which is later used for analysis and optimization. While we do provide the reader with a mathematical background, it is not required from a user to go in-depth into the inner workings of the inverse Laplace transform.

3.2. Description of the Algorithm

Mathematical optimization is a process of finding the best selection from a set of available options such as minimizing a function by choosing different input values. Usually, the input values are bounded. The type of domains and criteria for the best choice can vary largely depending on the type of problem. Therefore, these tasks are a significant part of applied mathematics. The algorithm starts with a calculation of random controller parameters within the limits.

• Let M be the number of simulations we would like to perform with a given transfer function.

• Assume $(K_{p,min},K_{p,max})$, $(K_{p,min},K_{p,\max })$ and $(K_{p,min},K_{p,max})$ to be the limits of the controller parameters, then

  $K_p=K_{p,min}+(K_{p,max}-K_{p,min})·{\rho }_1,$

  (22)$K_i=K_{i,min}+(K_{i,max}-K_{i,min})·{\rho }_2,$

  $K_d=K_{d,min}+(K_{d,max}-K_{d,min})·{\rho }_3,$

  where ${\rho }_1,{\rho }_2$ and ${\rho }_3$ are uniformly distributed random numbers. This algorithm can scan any region of parameter space of interest, however for the initial search it is safe to assume that $K_{p,min}=K_{i,min}=K_{d,min}=0.$ The upper limit, however, requires more attention since the equipment at hand may not handle too large values of the controller parameters. For example, in our case the upper limit of the electrical current in the coil was around 1A, or, in terms of voltage—$16.3\mathrm{V}$. It is easy to calculate the upper limits with the inverse Laplace transform of the expression

  (23)$\frac{1}{\Delta zs}\left(K_p+\frac{K_i}{s}+K_{d}s\right)\frac{1}{Ls+R},$

  where $\Delta z$ is the height (difference) a magnet needs to move up by. By substituting the $t=0$ in the resulting expression of $f\left(t\right)$, one may find whether the modeled response would correspond to that of a real plant.

For the given set of controller parameters, the algorithm continues by finding the roots ${\alpha }_m$ of the polynomial in the denominator $P\left(s\right)$, which in our case is $s(a_3{s}^4+a_2{s}^3+(a_1+A{K}_d)s^2+(a_0+A{K}_p)s+A{K}_i).$ If any of the roots are in the right-half plane, then the solution is unstable and we move to a new simulation. However, if all the roots are in the right-half plane the algorithm proceeds with the analysis of the resulting step-response signal. Before we start, let us introduce the variables and their initial values (see Figure 11).

• t– current time that the algorithm uses for calculations. Starts from zero;

• $f\left(t\right)$—the value of the step response signal at the current time t;

• $f^{\prime }\left(t\right)$—the value of the step response signal’s derivative at the current time t;

• $\Delta t$—the time step. This value changes as the algorithm progresses, depending on the step response. You can safely start with a very low value as the algorithm quickly finds the suitable time step for your process. For example, if one can ignore the changes occurring on the time scale of 1 ns, then the initial value can be $\Delta t=1$ ns.

• $\tilde{t}$—key parameter, the current recorded reference time. The two mechanisms for when this value changes are described later;

• $n=0$—how many decimal amplitude values the signal has crossed. An auxiliary counter needed for the first mechanism for detecting $\tilde{t}$;

• $k=0$—an auxiliary counter needed for the second mechanism for detecting $\tilde{t}$. Number of local extremums;

• $t_{min}$—the minimum out of all recorded reference times $\tilde{t}$. Needed for the calculation of the time step $\Delta t$. This parameter helps the algorithm to distinguish the possible rapid oscillations (see Figure 12);

• $N=100$—a positive integer regulating how fine do we divide the shortest reference time $t_{min}$ to calculate the $\Delta t$. Increase this parameter for additional accuracy;

• $t_{max}$—the largest reference time $\tilde{t}$ used as a verification for determining the settling time $T_s$;

• $t_k^{ex}$—the time of a local extremum;

• $f_k^{peak}$—the amplitude of a local extremum;

• $t_{3\%}$—the time when the signal reached the $\pm 3\%$ range;

• $T_s=0$—the settling time of the current curve;

• $PO$—percentage overshoot;

• $O{A}_k$—the amplitude difference between neighboring oscillations;

• $O{A}_{max}=0$—largest amplitude between neighboring oscillations. This helps to understand the scale of oscillations of the signal. If the number of local extremums $k\le 1,$ then this value is not defined.

Now, using pseudo code, let us describe the inner workings of the algorithm (see Algorithm 1). The algorithm scans the controller parameters’ space while performing the fast inverse Laplace transform calculations. It provides us with an explicit formula for the response of the initial system to a unit step signal. At the same time, it determines important signal parameters which are the number of oscillations before settling time, peak value, overshoot peak time, etc. For example, these signal data allow us to sort out highly oscillating responses.

One important feature of the designed algorithm is the fact that it automatically takes into account the Nyquist–Kotelnikov–Shannon theorem. This theorem specifies that a sinusoidal function in time or distance can be regenerated with no loss of information as long as it is sampled at a frequency greater than or equal to twice the frequency of oscillation. As shown above, the second mechanism of determining $\tilde{t}$ ensures the appropriate level of time resolution. In is important to keep this in mind while using some of the commercially available software. Let us show how, using a software without directly controlling the sampling rate may result in a wrong representation of a step response of a given system. To demonstrate the announced effect, we crafted a special transfer function of the same type as before (21), but $A=4.028\times 4.86\times {10}^5,a_3=1,a_2=207.7,a_1=-1257,a_0=-2.61\times {10}^5;$ and $K_p=14,K_i=1.6,K_d=30.$ Assume one needs to know the step response characteristics of the current process represented as this transfer function.

For this demonstration, the time step $\Delta t$ was manually controlled (see Figure 12). If we select the time step to be larger, like on the first plot—the transient response is a smooth curve with little to no overshoot. However, once we start decreasing the time step, oscillations emerge. At first these oscillations are angular. This indicates that the time moments at which we calculate the signal response miss some of the oscillations. This becomes obvious when, after the time step of 81 microseconds (the last two plots), new oscillations stop emerging and no new peaks appear. This means that we have reached the needed accuracy to truly represent this system’s step response. On its own, the algorithm chose the appropriate time step in less than 81 microseconds (for this process the last recorded time step was $\Delta t=2$ microseconds). This feature makes the developed algorithm adaptive to different time scales.

Algorithm 1 Processing of the given signal. Using the following algorithm, we gather large statistical data

Part 1. while$T_s=0$▹ continue until the settling time $T_s$ is determined $t\leftarrow t+\Delta t.$▹ increasing time by $\Delta t$ The first mechanism of recording $\tilde{t}$ (orange dots on Figure 11): if$n\le \left[f\left(t\right)/0.1\right]<n+1$then▹ where [ ] represents the floor function     $n++;\tilde{t}=(2t+\Delta t)/2;$ ▹ when the signal crosses decimal values a reference time is recorded (orange dots on Figure 11).     if $n+k=1$ then ▹ if this is the first a reference time $\tilde{t}$ is recorded, then         $t_{min}=\tilde{t};t_{max}=\tilde{t};$ ▹ assign this value to both $t_{min}$ and $t_{max}$     else         $t_{max}=\tilde{t};$ ▹ renewing only the maximum reference time     end if end if $\Delta t=t_{min}/N;$▹ calculate time step See Part 2 for further explanation. Part 2. The second mechanism of recording $\tilde{t}$ involves an extremum search (green dots on Figure 11). Since the function $f\left(t\right)$ is continuous, a moment of time when the derivative $f^{\prime }\left(t\right)$ changes sign implies a local extremum. if$f^{\prime }\left(t\right)·f^{\prime }(t-\Delta t)<0$then▹$f^{\prime }\left(t\right)$ changes sign, so a local extremum occurs in $(t-\Delta t,t)$     $t_k^{ex}=(2t+\Delta t)/2;k++;$ ▹ record the local extremum time and the number of them     $f_k^{peak}=f\left(t_k^{ex}\right);$ ▹ amplitude at the extremum     if $f_k^{peak}\ge f_{k-1}^{peak}$ then ▹ we have found a larger peak         $PO=(f_k^{peak}-1)·100\%$ ▹ percentage overshoot     end if     if $k>1$ then ▹ if there are multiple extremums         $O{A}_k=|f_k^{peak}-f_{k-1}^{peak}|$ ▹ the amplitude difference between the neighbor oscillations         $O{A}_{max}=max(O{A}_k,O{A}_{max})$ ▹ recording the maximum $OA$     end if     if $k=1$ then ▹ if this is the first extremum. This is done so that the $t_{max}$ does not increase uncontrollably.         $t_{max}=max(t_{max},t_k^{ex});$ ▹ the condition $k=1$ is needed so that $t_{max}$ does not increase uncontrollably     else         $\tilde{t}=t_k^{ex}-t_{k-1}^{ex};$▹ record the time between oscillations for additional accuracy (Nyquist theorem)         $t_{min}=min(t_{min},\tilde{t});t_{max}=max(t_{max},\tilde{t});$ ▹ in case this time is lower that the current $t_{min}$ or larger than the current value of $t_{max}$     end if end if $\Delta t=t_{min}/N;$▹ calculate time step Cont. At last, let us describe the test for finding $T_s.$ if$f\left(t\right)\in [0.97,1.03]$then▹ the signal appeared in the 3% range     $t_{3\%}=t;m=0;$ ▹ record the time and assign zero to the logic parameter m needed for the following cycle     while $t\le t_{3\%}+t_{max}$ AND $m=0$ do ▹ starting from $t_{3\%}$ for the time length of $t_{max}$ we check for the 3% criteria         $t=t+\Delta t$         if $f\left(t\right)\notin [0.97,1.03]$ then ▹ the signal exited the 3% area            $m=-1;$ ▹ therefore, the test for settling time stops and the algorithm resumes its work         end if     end while     if $m\ne -1$ then ▹ if the 3% condition did not break once during this test         $T_s=t_{3\%}$ ▹ we have found the settling time     end if end if

4. Results and Discussion

4.1. Multi-Objective Optimization

The developed algorithm returns the values of settling time $T_s$, percentage overshoot $PO$, number of local extremums k and the largest amplitude between the oscillations $O{A}_{max}$ in the form of Table 4. Using this table, one can easily arrange the column with a certain step-response characteristic (overshoot, for example) in ascending order to chose those controller parameters that yielded the preferred value (in our case, the lowest overshoot). The user is provided with enough information to sort out oscillating responses or understand their magnitude, using the parameter $O{A}_{max}$. The reason we try to avoid too much overshoot and oscillations is due to the nonlinearity of the system. A large overshoot may result in the inadequacy of the linearized model and instability of the real system.

Essentially, most controller tuning or designing processes are multi-objective problems. In our case, the first objective function is the overshoot $PO$ and the second one is the settling time $T_s$ while the PID parameters $(K_p,K_i,K_d)$ are selected as problem decision variables. The properties of magnetic levitation depend on parameters of the controller. In industry settings, controller parameters obtained using multi-objective genetic algorithms lead to a reduction of energy consumption, improved performance overall or and increased production rate [51,52].

In order to compare the effectiveness of the developed algorithm, we used common minimum search methods. One of them is the NSGA-II (Non-dominated Sorting Genetic Algorithm) with overshoot $PO,\%$ and settling time $T_s$ as the two objective functions. NSGA-II deals with a set of solutions simultaneously which improves the computational speed. Quite often this feature allows it to select several solutions of Pareto set in a single run of the algorithm [53].

To investigate this further, a search using the PSO (Particle Swarm Optimization) algorithm was performed with overshoot $PO,\%$ being the single objective function. PSO is an optimization method, which iteratively tries to improve a solution with regard to a given measure of quality. A particle is usually an element in some vector space—in our case $(K_p,K_i,K_d)$. PSO performs searching via a swarm of particles that updates each iteration. Using simple formulae, each particle moves in a direction depending on its previous best position and the best position among all of the particles in the whole swarm. Optimization ends if the relative change in the objective value over the last iterations is less than a defined function tolerance [54].

For the NSGA-II and PSO the Global Optimisation Toolbox in Matlab has been used. The decision variables of the search (the controller parameters) were limited to $K_p\in [0,200],K_p\in [0,250],K_p\in [0,10]$ due to hardware limitations (23). Figure 13 shows that NSGA-II outputs one point on the Pareto-front. This result stays the same after around one thousand total function evaluations. Usually in these situations, a rule of thumb is to find the single-objective minima to start the algorithm search from. By doing so, one may be able to obtain a wider Pareto set by avoiding a potential local minima. It is possible, however, that one still ends up with a single point, which means there is only one feasible Pareto point. In this study, starting from various points in the controller parameter space did not yield a different result.

The search using PSO yielded the same solution as seen in Figure 13. With this, we may conclude that there is no tradeoff curve (Pareto front) because all the objectives are minimized at the same point. The developed algorithm, on the other hand, was not affected by any local minima. While it is not a genetic-based algorithm the computational time was around the same as for the other two methods (simulations showed that all of the three algorithms required around ${10}^3$ function evaluations).

4.2. Experimental Verification

To provide an experimental confirmation, as well as to prove the effectiveness of the designed algorithm, we compared its results to the measurements using the magnetic levitation system, thoroughly described in the corresponding section. First, we made the magnet levitate steadily at $z_0=-25$ mm with an empirical tune of $(K_p,K_i,K_d)=(150,45,6.25)$. There are two main reasons behind this. First is to standardize the initial conditions before the start of each experiment. Second is for the system to be more comparable with the linearization. Then, we inserted new PID parameters in the microcontroller software to make sure the system is still stable. If it was, we forced the permanent magnet to rise $\Delta z=+1$ mm, while recording the transient response.

The comparison between the experiments and the simulations performed by the algorithm can be seen in Figure 14. The transients are between $-25$ mm and $-24$ mm. Both states are very close to the point of linearization, so the system performs fairly well and the results agree with the modeling. The relative values of overshoots related to the measured curves correspond with those obtained by the algorithm. The step responses with larger overshoot are expected to diverge more from the modeled ones. The divergence increases as the permanent magnet travels further from the equilibrium point at which the magnetic force function was expanded into the Taylor series.

For this research, as was shown before, Formula (5) has an approximation with the sum of squared errors lower by several orders of magnitude than the most commonly used formulas. However, the linearization still adds some inaccuracy. The real magnetic force changes with distance z by a different law. It is quite difficult to carefully estimate a finite size coil’s magnetic force applied to a permanent magnet, when the dimensions of the coil are comparable to the distance to the permanent magnet. Another obstacle would be taking into account the shape of the permanent magnet. These are serious mathematical problems that go beyond this research’s goal which is to demonstrate the capabilities of the designed controller tuning and system design algorithm. For the purpose of this research, it is sufficient to use the approximation we presented in Section The Magnetic Levitation System—Figure 14 reflects that. We also discuss the possibility of improving the accuracy of the modeling in the following section.

4.3. The Multi-Objective Problem for Optimal Coil Parameters

Here we show how, by taking into account more information about the magnetic levitation system and varying one or more of the coil’s parameters, one may find their optimal values. It is possible to use the data provided by the algorithm as a criterion for optimization. One may look at it as a Monte-Carlo method of sorts. Here we discuss an example of this method’s realization. Before, as an approximation we used formula (5) for the magnetic field of the coil:

$B_z=\frac{{\mu }_0}{2\pi }i\frac{\pi R^2}{{(z^2+R^2)}^{3/2}}.$

This is, however, the magnetic field caused by just one of the rings within the coil. For that reason, let us perform an integration over the length of the coil. We should point out, though, that any mistake in the formulas at this stage will result in a wrong analysis overall. So extra caution should be taken to keep the expressions correct. In our setup, the permanent magnet is situated below the coil, which means $z<0$. That is why

$B_{zcoil}=\frac{{\mu }_0}{2\pi }ni{\int }_0^L\frac{\pi R^{2}dl}{{({(l-z)}^2+R^2)}^{3/2}}=$

(24)$=\frac{{\mu }_0}{2}ni\left(\frac{L-z}{{(L-z)}^2+R^2}+\frac{z}{z^2+R^2}\right),$

where L is the length of the coil and n is the density of winding turns per unit length. To get the explicit expression of the magnetic force in (4), the final step is to perform the differentiation of this expression.

(25)${\frac{\partial B_z}{\partial z}}_{coil}=\frac{{\mu }_0}{2}ni\left(\frac{{(L-z)}^2-R^2}{{({(L-z)}^2+R^2)}^2}+\frac{R^2-z^2}{{(z^2+R^2)}^2}\right),$

The linearization of the process involved requires a Taylor transformation, after which Newton’s equation is converted to

(26)$m\frac{d^2\widehat{z}}{d{t}^2}=\frac{\partial F_z}{\partial z}(i_0,z_0)·\widehat{z}+\frac{\partial F_z}{\partial i}(i_0,z_0)·\widehat{i},$

where

$\frac{\partial F_z}{\partial z}(i_0,z_0)=-2{\mu }_z\frac{{\mu }_0}{2}n{i}_0\times$

$\times \left(\frac{(L-z_0)(3R^2-{(L-z_0)}^2)}{{({(L-z_0)}^2+R^2)}^3}+\frac{z_0(3R^2-z_0^2)}{{(z_0^2+R^2)}^3}\right),$

$\frac{\partial F_z}{\partial i}(i_0,z_0)={\mu }_z\frac{{\mu }_0}{2}n\times$

(27)$\left(\frac{{(L-z_0)}^2-R^2}{{({(L-z_0)}^2+R^2)}^2}+\frac{R^2-z_0^2}{{(z_0^2+R^2)}^2}\right).$

Now the acquired transfer function has the length of the coil L, the radius of the coil R and the density of winding turns per unit length n within the transfer function. This opens a whole new mathematical problem of optimizing these parameters for best performance. For this, our algorithm with its rapid calculations is a great tool for optimization. This simulation data can be the basis for a multi-objective optimization to determine the best values for these parameters. First, we postulate an optimization criterion which may consist of one feature of the step response or several of them. Then, the algorithm provides characteristics of thousands of modeled step responses that can be visually represented as in Figure 13 using radar charts or a Pareto-optimal front.

5. Conclusions

In this study, we presented a framework for designing control systems. For this purpose, an algorithm was created which determines the key features of simulated transient responses. With modern computing power, this algorithm is suitable for multi-objective optimization tasks such as nominal parameters of complicated mechanical systems. The expected curve form of regular stable transient responses opens a possibility of fast curve analysis using the developed algorithm. This, in turn, opens the possibility of collecting large simulation data samples, which can be applied to solve different multi-optimization problems such as optimal coil parameters. A comparison with common genetic search methods showed the competitiveness of the designed algorithm.

With the crafted magnetic levitation system as a means of testing, we were able to show the algorithm’s versatility. Specific features were discussed, such as linearizing the magnetic levitation force on the permanent magnet as well as curve fitting of the measured data. In this case, the main focus was reducing the overshoot, for which the proposed method was able to provide numerous suitable sets of controller parameters. The transient response behavior, as well as values of overshoot related to the measured curves, corresponded well with those obtained by the simulation. Therefore, this method proved to be highly effective for pre-production purposes.

The main benefits of the developed algorithm can be summarized. Coding, applicable for any programming language. An important feature is its immunity to the effect of local minima trapping. As shown in Figure 13, the developed algorithm in the given conditions proved to be versatile and provided multiple suitable solutions to the multi-objective problem between objective functions of overshoot and settling time. In addition, we showed how the algorithm automatically adjusts the size of the time step of the simulated signal to meet the conditions of the Nyquist–Kotelnikov–Shannon theorem. As we mentioned, this feature is particularly important for open-loop unstable systems (such as the magnetic levitation system) since unsupervised large oscillations may result in unwanted sideways irregularities in motion and instability of the whole controlled system. This method is a basis for solving a multi-objective problem of optimal coil dimensions and proportions. Additional accuracy in the nonlinearity identification process will lead to a wider optimization problem of not only the controller parameters but also key features of any magnetic levitation system.

Author Contributions

Conceptualization, I.R. and P.P.; methodology, I.R. and P.P.; software, I.R.; validation, I.R. and P.P.; formal analysis, I.R.; investigation, I.R.; resources, I.R. and P.P.; data curation, I.R. and P.P.; writing—original draft preparation, I.R. and P.P.; writing—review and editing, I.R.; visualization, I.R. and P.P.; supervision, P.P. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used in this research is available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Two possible levitation setups (attraction-(left), repulsion-(right)) [9,14].

Figure 2. Scheme of a magnetic levitation system.

Figure 3. The basic block diagram of the system.

Figure 4. Experimental determination of the magnetic force.

Figure 5. Input voltage [Forumla omitted. See PDF.] and output voltage [Forumla omitted. See PDF.] in a scheme to determine the coil’s inductance [Forumla omitted. See PDF.].

Figure 6. Input voltage [Forumla omitted. See PDF.] (black) and output voltage [Forumla omitted. See PDF.] (blue) for the circuit shown in Figure 5.

Figure 7. The scheme of the actuator.

Figure 8. Possible sensors for detecting position of the permanent magnet.

Figure 9. Schematic representation of a coil.

Figure 10. The role of the microcontroller (Arduino) in the system.

Figure 11. The algorithm analyses a simulated transient response. The orange points represent the times when the output signal crosses a new reference value. The green points are the signal’s extrema. The key parameter of this process is the current recorded reference time [Forumla omitted. See PDF.]. The two mechanisms for when this value changes are described (see the algorithm’s description).

Figure 12. A showcase of modelling a mechanical system without satisfying the Nyquist theorem. For some mechanical systems it may be of great importance to distinguish oscillatory step responses. Such oscillations may, for example, result in unnecessary sideways oscillations which may result in complete loss of control over the object. The developed algorithm, however, requires no prior knowledge about the time scale of the process involved as it determines the appropriate time step [Forumla omitted. See PDF.] microseconds automatically.

Figure 13. A comparison between the developed algorithm and genetic search methods. With the used boundaries for the decision variables [Forumla omitted. See PDF.] the developed algorithm proved to be insensible to the local minima trapping effect as it showed multiple solutions on the Pareto set. Combined with fast computations, this makes for an optimal choice for a design method involving transient response analysis.

Figure 14. A comparison of a sequence of measured and modeled responses with reducing overshoot. Overall, the modeling is consistent with the experiments. The key feature influencing the compliance of the data is the accuracy of the used expression for the magnetic force. As expected, with larger values of overshoot (blue and yellow curves) the model becomes less accurate. The difficulties involving the nonlinearity identification of the magnetic force of the coil applied to the permanent magnet in a form of a small cylinder are discussed in the corresponding section.

Table 1

Magnetic force $F_m$ in millinewtons $\left[\mathrm{mN}\right]$ in relationship to the current i in Amperes $\left[\mathrm{A}\right]$ and the permanent magnet position z in millimeters $\left[\mathrm{mm}\right]$.

i A	0 A	0.2 A	0.4 A	0.6 A	0.8 A	1.0 A
$z=-15$ mm	0 mN	52.0 mN	101.0 mN	151.1 mN	200.1 mN	250.2 mN
$z=-20$ mm	0 mN	32.4 mN	63.8 mN	94.2 mN	126.5 mN	156.0 mN
$z=-25$ mm	0 mN	23.5 mN	45.1 mN	65.7 mN	86.3 mN	106.9 mN
$z=-30$ mm	0 mN	14.7 mN	28.4 mN	42.2 mN	54.9 mN	68.7 mN
$z=-35$ mm	0 mN	9.8 mN	18.6 mN	28.4 mN	38.3 mN	48.1 mN
$z=-40$ mm	0 mN	6.9 mN	13.7 mN	20.6 mN	27.5 mN	33.4 mN

Table 2

The results of data analysis using sum of squared errors.

Function	$-\mathit{aiz}/{({\mathit{z}}^2+{\mathit{R}}^2)}^{5/2}$	$\mathit{ai}/{\mathit{z}}^2$	$\mathit{ai}/\mathit{z}$	$\mathit{ai}/{\mathit{z}}^3$
SSE	$5.6\times {10}^{-5}$	0.0026	0.017	0.022

Table 3

Sensor circuit output $U_s\left[\mathrm{V}\right]$ in relationship to the permanent magnet position $z\left[\mathrm{mm}\right]$.

z [mm]	${\mathit{U}}_{\mathit{s}}$ [V]
$-15$	3.32
$-20$	1.20
$-25$	0.504
$-30$	0.219
$-35$	0.135
$-40$	0.072

Table 4

Step-response characteristics provided by the tuning algorithm. The value of $O{A}_{max}$ provides us with the scale of oscillations. If the number of local extremums $k\le 1,$ then this value is obviously not defined.

${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{d}}$	${\mathit{T}}_{\mathit{s}},$ (s)	k	${\mathit{OA}}_{\mathit{max}}$	PO	${\mathit{t}}_{\mathit{min}},$ (s)
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
223	875	9.58	0.49	1	—	61.3	0.0026
534	938	6.75	0.74	3	0.59	80.7	0.0032
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
85.1	171	3.67	1.41	4	2.7	129	0.0047
224	391	8.15	0.26	1	—	63.5	0.0030
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

References

1. Astrom, K.J.; Hagglund, T. Advanced PID Control; ISA-The Instrumentation, Systems, and Automation Society: Research Triangle Park, NC, USA, 2006.

2. Ang, K.H.; Chong, G.; Li, Y. PID control system analysis, design, and technology. IEEE Trans. Control. Syst. Technol.; 2005; 13, pp. 559-576. [DOI: https://dx.doi.org/10.1109/TCST.2005.847331]

3. Astolfi, A. Model reduction by moment matching for nonlinear systems. Proceedings of the 2008 47th IEEE Conference on Decision and Control; Cancún, Mexico, 9–11 December 2008; pp. 4873-4878. [DOI: https://dx.doi.org/10.1109/CDC.2008.4738791]

4. Shrivastava, N.; Varshney, P. Comparative analysis of order reduction techniques. Proceedings of the 2016 Second International Innovative Applications of Computational Intelligence on Power, Energy and Controls with their Impact on Humanity (CIPECH); Ghaziabad, India, 18–19 November 2016; pp. 46-50. [DOI: https://dx.doi.org/10.1109/CIPECH.2016.7918735]

5. Doetsch, G. Introduction to the Theory and Application of the Laplace Transformation; Springer: Berlin/Heidelberg, Germany, 1974.

6. Cohen, A. Numerical Methods for Laplace Transform Inversion; Numerical Methods and Algorithms; Springer: Berlin/Heidelberg, Germany, 2007.

7. Davies, B.; Martin, B. Numerical inversion of the Laplace transform: A survey and comparison of methods. J. Comput. Phys.; 1979; 33, pp. 1-32. [DOI: https://dx.doi.org/10.1016/0021-9991(79)90025-1]

8. Hosono, T. Numerical inversion of Laplace transform and some applications to wave optics. Radio Sci.; 1981; 16, pp. 1015-1019. [DOI: https://dx.doi.org/10.1029/RS016i006p01015]

9. Yan, L. Development and application of the maglev transportation system. IEEE Trans. Appl. Supercond.; 2008; 18, pp. 92-99.

10. Wai, R.J.; Lee, J.D.; Chuang, K.L. Real-time PID control strategy for maglev transportation system via particle swarm optimization. IEEE Trans. Ind. Electron.; 2010; 58, pp. 629-646. [DOI: https://dx.doi.org/10.1109/TIE.2010.2046004]

11. Peijnenburg, A.; Vermeulen, J.; Van Eijk, J. Magnetic levitation systems compared to conventional bearing systems. Microelectron. Eng.; 2006; 83, pp. 1372-1375. [DOI: https://dx.doi.org/10.1016/j.mee.2006.01.248]

12. Fang, J.; Le, Y.; Sun, J.; Wang, K. Analysis and design of passive magnetic bearing and damping system for high-speed compressor. IEEE Trans. Magn.; 2012; 48, pp. 2528-2537. [DOI: https://dx.doi.org/10.1109/TMAG.2012.2196443]

13. Barry, N.; Casey, R. Elihu Thomson’s jumping ring in a levitated closed-loop control experiment. IEEE Trans. Educ.; 1999; 42, pp. 72-80. [DOI: https://dx.doi.org/10.1109/13.746338]

14. Boudali, H.; Williams, R.; Giras, T. A Simulink simulation framework of a MagLev model. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit; 2003; 217, pp. 227-236. [DOI: https://dx.doi.org/10.1243/095440903769012911]

15. Berkelman, P.; Dzadovsky, M. Magnetic levitation over large translation and rotation ranges in all directions. IEEE/ASME Trans. Mechatronics; 2011; 18, pp. 44-52. [DOI: https://dx.doi.org/10.1109/TMECH.2011.2161614]

16. Yaseen, M.H.; Abd, H.J. Modeling and control for a magnetic levitation system based on SIMLAB platform in real time. Results Phys.; 2018; 8, pp. 153-159. [DOI: https://dx.doi.org/10.1016/j.rinp.2017.11.026]

17. Chopade, A.S.; Khubalkar, S.W.; Junghare, A.; Aware, M.; Das, S. Design and implementation of digital fractional order PID controller using optimal pole-zero approximation method for magnetic levitation system. IEEE/CAA J. Autom. Sin.; 2016; 5, pp. 977-989. [DOI: https://dx.doi.org/10.1109/JAS.2016.7510181]

18. Morales, R.; Sira-Ramírez, H. Trajectory tracking for the magnetic ball levitation system via exact feedforward linearisation and GPI control. Int. J. Control; 2010; 83, pp. 1155-1166. [DOI: https://dx.doi.org/10.1080/00207171003642196]

19. Awelewa, A.A.; Samuel, I.A.; Abdulkareem, A.; Iyiola, S.O. An Undergraduate Control Tutorial on Root Locus-Based Magnetic Levitation System Stabilization. Int. J. Eng. Comput. Sci.; 2013; 13, pp. 22-30.

20. Hurley, W.G.; Wolfle, W.H. Electromagnetic design of a magnetic suspension system. IEEE Trans. Educ.; 1997; 40, pp. 124-130. [DOI: https://dx.doi.org/10.1109/13.572325]

21. Oguchi, K.; Tomigashi, Y. Digital control for a magnetic suspension system as an undergraduate project. Int. J. Electr. Eng. Educ.; 1990; 27, pp. 226-236. [DOI: https://dx.doi.org/10.1177/002072099002700305]

22. Wei, W.; Xue, W.; Li, D. On disturbance rejection in magnetic levitation. Control Eng. Pract.; 2019; 82, pp. 24-35. [DOI: https://dx.doi.org/10.1016/j.conengprac.2018.09.018]

23. Salim, T.T.; Karsli, V.M. Control of single axis magnetic levitation system using fuzzy logic control. Int. J. Adv. Comput. Sci. Appl.; 2013; 4, pp. 83-88.

24. Al-Muthairi, N.; Zribi, M. Sliding mode control of a magnetic levitation system. Math. Probl. Eng.; 2004; 2004, pp. 93-107. [DOI: https://dx.doi.org/10.1155/S1024123X04310033]

25. Kuo, C.L.; Li, T.H.S.; Guo, N.R. Design of a novel fuzzy sliding-mode control for magnetic ball levitation system. J. Intell. Robot. Syst.; 2005; 42, pp. 295-316. [DOI: https://dx.doi.org/10.1007/s10846-004-3026-3]

26. Oliveira, V.A.; Tognetti, E.S.; Siqueira, D. Robust controllers enhanced with design and implementation processes. IEEE Trans. Educ.; 2006; 49, pp. 370-382. [DOI: https://dx.doi.org/10.1109/TE.2006.879263]

27. de Jesús Rubio, J.; Zhang, L.; Lughofer, E.; Cruz, P.; Alsaedi, A.; Hayat, T. Modeling and control with neural networks for a magnetic levitation system. Neurocomputing; 2017; 227, pp. 113-121. [DOI: https://dx.doi.org/10.1016/j.neucom.2016.09.101]

28. Hajjaji, A.E.; Ouladsine, M. Modeling and nonlinear control of magnetic levitation systems. IEEE Trans. Ind. Electron.; 2001; 48, pp. 831-838. [DOI: https://dx.doi.org/10.1109/41.937416]

29. Morales, R.; Feliu, V.; Sira-Ramirez, H. Nonlinear control for magnetic levitation systems based on fast online algebraic identification of the input gain. IEEE Trans. Control Syst. Technol.; 2011; 19, pp. 757-771. [DOI: https://dx.doi.org/10.1109/TCST.2010.2057511]

30. Truong, T.N.; Vo, A.T.; Kang, H.J. Real-Time Implementation of the Prescribed Performance Tracking Control for Magnetic Levitation Systems. Sensors; 2022; 22, 9132. [DOI: https://dx.doi.org/10.3390/s22239132] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36501837]

31. Maximov, S.; Gonzalez-Montañez, F.; Escarela-Perez, R.; Olivares-Galvan, J.C.; Ascencion-Mestiza, H. Analytical Analysis of Magnetic Levitation Systems with Harmonic Voltage Input. Actuators; 2020; 9, 82. [DOI: https://dx.doi.org/10.3390/act9030082]

32. Tipler, P.A.; Mosca, G. Physics for Scientists and Engineers; 6th ed. W. H. Freeman and Company: New York, NY, USA, 2008.

33. Serway, R.A.; Jewett, J.W. Physics for Scientists and Engineers with Modern Physics; 9th ed. Brooks/Cole: Salt Lake City, UT, USA, 2014.

34. de Queiroz, A.C.M. Mutual Inductance and Inductance Calculations by Maxwell’s Method. 2014; Available online: http://www.coe.ufrj.br/~acmq/programs (accessed on 7 March 2020).

35. Galvão, R.K.H.; Yoneyama, T.; Araújo, F.M.U.; Machado, R.G. A simple technique for identifying a linearized model for a didactic magnetic levitation system. IEEE Trans. Educ.; 2003; 46, pp. 22-25. [DOI: https://dx.doi.org/10.1109/TE.2002.804403]

36. Lin, C.M.; Lin, M.H.; Chen, C.W. SoPC-based adaptive PID control system design for magnetic levitation system. IEEE Syst. J.; 2011; 5, pp. 278-287. [DOI: https://dx.doi.org/10.1109/JSYST.2011.2134530]

37. Zhang, Y.; Xian, B.; Ma, S. Continuous robust tracking control for magnetic levitation system with unidirectional input constraint. IEEE Trans. Ind. Electron.; 2015; 62, pp. 5971-5980. [DOI: https://dx.doi.org/10.1109/TIE.2015.2434791]

38. Bächle, T.; Hentzelt, S.; Graichen, K. Nonlinear model predictive control of a magnetic levitation system. Control Eng. Pract.; 2013; 21, pp. 1250-1258. [DOI: https://dx.doi.org/10.1016/j.conengprac.2013.04.009]

39. Folea, S.; Muresan, C.I.; De Keyser, R.; Ionescu, C.M. Theoretical analysis and experimental validation of a simplified fractional order controller for a magnetic levitation system. IEEE Trans. Control Syst. Technol.; 2015; 24, pp. 756-763. [DOI: https://dx.doi.org/10.1109/TCST.2015.2446496]

40. Lundberg, K.H.; Lilienkamp, K.A.; Marsden, G. Low-cost magnetic levitation project kits. IEEE Control Syst. Mag.; 2004; 24, pp. 65-69.

41. Chalupa, P.; Maly, M.; Novák, J. Nonlinear Simulink Model Of Magnetic Levitation Laboratory Plant. Proceedings of the ECMS; Regensburg, Germany, 31 May–3 June 2016; pp. 293-299.

42. Chalupa, P.; Novák, J.; Malỳ, M. Modelling and model predictive control of magnetic levitation laboratory plant. Proceedings of the 31st European Conference on Modelling and Simulation, ECMS 2017; Budapest, Hungary, 23–26 May 2017; European Council for Modelling and Simulation: Caserta, Italy, 2017.

43. da Silva, L.R.; Flesch, R.C.C.; Normey-Rico, J.E. Controlling industrial dead-time systems: When to use a PID or an advanced controller. ISA Trans.; 2020; 99, pp. 339-350. [DOI: https://dx.doi.org/10.1016/j.isatra.2019.09.008]

44. Sánchez, J.; Guinaldo, M.; Visioli, A.; Dormido, S. Identification of process transfer function parameters in event-based PI control loops. ISA Trans.; 2018; 75, pp. 157-171. [DOI: https://dx.doi.org/10.1016/j.isatra.2018.01.033] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29455888]

45. Balaguer, P.; Alfaro, V.; Arrieta, O. Second order inverse response process identification from transient step response. ISA Trans.; 2011; 50, pp. 231-238. [DOI: https://dx.doi.org/10.1016/j.isatra.2010.11.005]

46. Ozsoy, M.; Sims, N.D.; Ozturk, E. Robotically assisted active vibration control in milling: A feasibility study. Mech. Syst. Signal Process.; 2022; 177, 109152. [DOI: https://dx.doi.org/10.1016/j.ymssp.2022.109152]

47. Chen, J.; Ma, D.; Xu, Y.; Chen, J. Delay Robustness of PID Control of Second-Order Systems: Pseudoconcavity, Exact Delay Margin, and Performance Tradeoff. IEEE Trans. Autom. Control; 2022; 67, pp. 1194-1209. [DOI: https://dx.doi.org/10.1109/TAC.2021.3059155]

48. Mercorelli, P.; Lehmann, K.; Liu, S. Robust flatness based control of an electromagnetic linear actuator using adaptive PID controller. Proceedings of the 42nd IEEE International Conference on Decision and Control (IEEE Cat. No.03CH37475); Maui, HI, USA, 9–12 December 2003; Volume 4, pp. 3790-3795. [DOI: https://dx.doi.org/10.1109/CDC.2003.1271739]

49. Mercorelli, P. An Antisaturating Adaptive Preaction and a Slide Surface to Achieve Soft Landing Control for Electromagnetic Actuators. IEEE/ASME Trans. Mechatronics; 2012; 17, pp. 76-85. [DOI: https://dx.doi.org/10.1109/TMECH.2010.2089467]

50. Zhao, C.; Guo, L. Control of Nonlinear Uncertain Systems by Extended PID. IEEE Trans. Autom. Control; 2021; 66, pp. 3840-3847. [DOI: https://dx.doi.org/10.1109/TAC.2020.3030876]

51. Babajamali, Z.; Khabaz, M.K.; Aghadavoudi, F.; Farhatnia, F.; Eftekhari, S.A.; Toghraie, D. Pareto multi-objective optimization of tandem cold rolling settings for reductions and inter stand tensions using NSGA-II. ISA Trans.; 2022; 130, pp. 399-408. [DOI: https://dx.doi.org/10.1016/j.isatra.2022.04.002]

52. Yusoff, Y.; Ngadiman, M.S.; Zain, A.M. Overview of NSGA-II for Optimizing Machining Process Parameters. Procedia Eng.; 2011; 15, pp. 3978-3983. [DOI: https://dx.doi.org/10.1016/j.proeng.2011.08.745]

53. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput.; 2002; 6, pp. 182-197. [DOI: https://dx.doi.org/10.1109/4235.996017]

54. Mathworks. Global Optimization Toolbox User’s Guide; Mathworks: Natick, MA, USA, 2018.