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

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.


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].

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).  [9,14].
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.

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.

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 B (or its component B z ) of a single ring can be obtained by the application of the Biot-Savart law, where µ 0 = 4π · 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 µ = µ x , µ y , µ z and be in a uniform magnetic field B = B x , B y , B z . If the loop is small enough, then the torque acting on such a loop is given by a simple expression [33]: therefore, the force F acting on a magnetic dipole in this field is the gradient of the potential energy associated with this torque: As magnetic dipole moment µ has only a vertical component (µ z in this case), which is constant, we can write the vertical component of the force F z as follows: 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): where µ 0 = 4π · 10 −7 H / m is the magnetic permeability of vacuum, µ 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 mm, while the mean radius of the coil's winding is R = 37 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 |µ z | = 0.49 A · m 2 . As a means of comparison we chose the commonly used functions: 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 × 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.

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 m and the cross sectional area was 0.246 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 Ω. In order to determine the inductance of the coil L c , a resistor R r = 468 Ω was added in series to the coil as shown schematically in Figure 5. The ratio between output voltage U out (t) and input voltage U in (t) can be obtained by a simple voltage divider rule, The ratio of amplitudes of the output and input voltages is, therefore, given by: By a simple algebraic manipulation, we get the following result: The actual waveforms for both U in (t) and U out (t) are shown in Figure 6. For the signal U in = 5.12 V and ω = 6280 rad/s we get the U out = 2.87 V. Finally, using Equation (9) we get L c = 52.1 mH.

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(t).
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: where m = 3 grams is the mass of the magnet, g is the gravitational acceleration and the force F z is given by the expression (5). Letẑ(t) be the relative change of coordinate z from the initial state z 0 : Zero initial valueẑ(+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 mm. Using (5) to calculate the initial current from equation We expand the function F z in a Taylor series at the point (i 0 , z 0 ), whereî = 0,ẑ = 0. Since the operating point is chosen in a way that mg = F z (i 0 , z 0 ) we rewrite the Equation (11) as where and The Laplace Transform of the Equation (13) is 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: Combined with the expression (10) the resulting transfer function of our magnetic levitation system is

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 PW M (which is the output voltage of the controller shown in Figure 7). The maximum value of U PW M is equal to 3.33 V. The voltage U cc was selected to be 16.2 V, so that the maximum current through the electromagnet is 1 A. The transfer function of the actuator is, therefore, equal to:

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 photoreceiver [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: 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 (U s ) in the controller.

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 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 Figures 3 and 10). The control algorithm is the common PID algorithm. The control signal (output of the controller) is a pulse width modulated signal.

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 stepresponse 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 (s) = G a (s) · G o (s), and G ao (s) = 3723 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 (s) · G −1 s (s) = 1, we have a unity feedback control system. The transfer function of the PID controller in a parallel form is G c (s) = 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(s) ends up being 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 (t) and its derivative f (t) 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.

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 ) · ρ 1 , where ρ 1 , ρ 2 and ρ 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 V. It is easy to calculate the upper limits with the inverse Laplace transform of the expression where ∆z is the height (difference) a magnet needs to move up by. By substituting the t = 0 in the resulting expression of f (t), 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 α m of the polynomial in the denominator P(s), which in our case is s(a 3 s 4 + a 2 s 3 + (a 1 + AK d )s 2 + (a 0 + AK p )s + AK 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). The key parameter of this process is the current recorded reference timet. The two mechanisms for when this value changes are described (see the algorithm's description).
• tcurrent time that the algorithm uses for calculations. Starts from zero; • f (t)-the value of the step response signal at the current time t; • f (t)-the value of the step response signal's derivative at the current time t; • ∆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 ∆t = 1 ns. •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 detectingt; • k = 0-an auxiliary counter needed for the second mechanism for detectingt. Number of local extremums; • t min -the minimum out of all recorded reference timest. Needed for the calculation of the time step ∆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 ∆t. Increase this parameter for additional accuracy; 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 determiningt 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 × 4.86 × 10 5 , a 3 = 1, a 2 = 207.7, a 1 = −1257, a 0 = −2.61 × 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 ∆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 ∆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 ← t + ∆t. increasing time by ∆t The first mechanism of recordingt (orange dots on Figure 11): if n ≤ [ f (t)/0.1] < n + 1 then where [ ] represents the floor function n + +;t = (2t + ∆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 timet is recorded, then t min =t; t max =t; assign this value to both t min and t max else t max =t; renewing only the maximum reference time end if end if ∆t = t min /N; calculate time step See Part 2 for further explanation. Part 2. The second mechanism of recordingt involves an extremum search (green dots on Figure 11). Since the function f (t) is continuous, a moment of time when the derivative f (t) changes sign implies a local extremum. if f (t) · f (t − ∆t) < 0 then f (t) changes sign, so a local extremum occurs in (t − ∆t, t) t ex k = (2t + ∆t)/2; k + +; record the local extremum time and the number of them f OA max = max(OA k , OA 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 ex k ); the condition k = 1 is needed so that t max does not increase uncontrollably elset = t ex k − t ex k−1 ; record the time between oscillations for additional accuracy (Nyquist theorem) t min = min(t min ,t); t max = max(t max ,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 ∆t = t min /N; calculate time step At last, let us describe the test for finding T s .

Algorithm 1 Cont.
if f (t) ∈ [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 ≤ 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 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 = −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 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 ∆t = 2 microseconds automatically.

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 OA max in the form of Table 4. Using this table, one can easily arrange the column with a certain stepresponse 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 OA 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. Table 4.
Step-response characteristics provided by the tuning algorithm. The value of OA max provides us with the scale of oscillations. If the number of local extremums k ≤ 1, then this value is obviously not defined.
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 ∈ [0, 200], K p ∈ [0, 250], K p ∈ [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).

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 ∆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.

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: 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 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.
The linearization of the process involved requires a Taylor transformation, after which Newton's equation is converted to where 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.

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.