Fractional PI^λD Controller Design for a Magnetic Levitation System

Abstract

Currently, there are no formalized methods for tuning non-integer order controllers. This is due to the fact that implementing these systems requires using an approximation of the non-integer order terms. The Oustaloup approximation method of the $s^{\alpha }$ fractional derivative is intuitive and widely adopted in the design of fractional-order PI${}^{\lambda }$D controllers. It requires special considerations for real-time implementations as it is prone to numerical instability. In this paper, for design and tuning of fractional regulators, we propose two methods.The first method relies on Nyquist stability criterion and stability margins. We base the second on parametric optimization via Simulated Annealing of multiple performance indicators. We illustrate our methods with a case study of the PI${}^{\lambda }$D controller for the Magnetic Levitation System. We illustrate our methods’ efficiency with both simulations and experimental verification in both nominal and disturbed operation.

1. Introduction

The magnetic levitation system is a highly relevant area of control design. It is because these systems have many varied uses such as in minimal friction bearings, high-speed maglev passenger trains, levitation of wind tunnel models, vibration isolation of sensitive machinery, levitation of molten metal in induction furnaces, and the levitation of metal slabs during manufacture, see [1]. Research on this topic has led to the development of efficient magnetic bearings (see [2]) and flywheels for space satellite low earth orbits, hybrid electric vehicles (see [3]), and many stationary applications. An important direction of research in this area is also an active magnetic bearing, which was described in their works by among others—Srinivas (see [4]), Uzhegov (see [5]), and Kandil (see [6]).

A very important area of the research about the active magnetic levitation system is control structure. They are naturally unstable, so they require control. Many researchers considered this topic (see [1]) and Joo and Seo (see [7]) designed linearizing feedback control for the magnetic levitation system to PID controllers. The comparison of this approach with Takagi–Sugeno fuzzy control (see [8]). Baranowski and Piątek (see [9]) presented a cascade variant of the linearizing feedback. Bloch (see [10]) discussed real-time neural feed-forward control. Piątek (see [11]) designed and implemented quick linear control for active magnetic levitation systems based on FPGA circuits.

Research on the control of magnetic levitation using non-integer (fractional) systems was one of the first attempts by Piłat (see [12], where a non-integer order PD controller was considered. Tepljakov et al. (see [13,14]) described the problem of fractional-order PID controller design for a model of a magnetic levitation system. The latest research focuses on the digital implementation of non-integer controller for a real plant; this topic was considered by Chopade et al. [15], Rojas et al. [16] and Ananthababu et al. [17]. Pandey et al. (see [18,19]) proposed an anti-windup fractional PID controller for magnetic levitation. Roy et al. (see [20]) presented a comparative study between the Sliding Mode Controller and Fractional Order Sliding Mode Controller when applied to a magnetic levitation system. Swain et al. (see [21]) presented real-time implementation of fractional order PID controllers for a magnetic levitation.

Most of these results that are simulated based or experimental are hard to replicate or use not precisely for fractional approximations, which makes the fractional part debatable.

Contributions of this paper are:

• we propose two different methods of controller parameter selection for the active magnetic levitation system,
• and we provide experimental verification of both design methods.

This paper is organized as follows. Section 2 presents a brief introduction to the mathematical model of the magnetic levitation laboratory system. In the next section, the authors present a non-integer PID controller and proposed two methods of tuning these parameters. In Section 4, the authors show an example of using a proposed method for construction of a PI${}^{\lambda }$D controller for the laboratory magnetic levitation system. The next section describes the implementation of the controller in real-time systems. The results of the experiments are described in Section 6.

2. The Mathematical Model of the Magnetic Levitation Laboratory System

The magnetic levitation system (see Figure 1) consists of the electromagnet, the ferromagnetic sphere (which we will be calling ”ball”), the current driver, and the position measurement system consisting of the light source and the light sensor.

The mathematical description of the system will be based on Newton’s second law, where $x_1\left(t\right)$ is the ball position, measured as the distance from the electromagnet, $x_3\left(t\right)$ is the electromagnet coil current, and $x_2\left(t\right)$ is the ball velocity. Ball speed can by calculated from value of force generated by the electromagnet, m is the mass of the ball, and g is the gravitational acceleration. The velocity of ball is given by formula (see [1,23]):

$$\begin{array}{ccc}\hfill {\dot{x}}_2\left(t\right)& =& -\frac{a}{2mb}exp\left(-\frac{x_1\left(t\right)}{b}\right)x_3^2+g\hfill \end{array}$$

(1)

where a and b are positive constants. Parameters a and b were determined by analysis of a series of steady state points of the system with a closed stabilizing feedback loop. The exponential function was fitted into these points through a least squares minimization. For details, see [23].

The coil current in the system usually is influenced by many factors like changes in inductance, velocity, and others. However, our system includes a current driver, which has its own feedback loop. This solution is very popular (see [24]) because it leads to either lower order or simpler model structure. In an optimal situation, the driver should allow full current control; however, in real situations, it introduces its own dynamics. For the system considered, these dynamics can be sufficiently modeled by a first order dynamical system given by the following equation:

$$\begin{array}{ccc}\hfill {\dot{x}}_3\left(t\right)& =& \frac{1}{T_s}(k_{s}u\left(t\right)-i_s-x_3\left(t\right))\hfill \end{array}$$

(2)

where $u\left(t\right)$ is the control voltage, $k_s$ is the gain of current controller, $T_s$ is the time constant of the current driver, and $i_s$ is the zero error of current driver.

Let us introduce state space vector $\mathbf{x}$ given by

$$\begin{array}{ccc}\hfill \mathbf{x}& =& {\left[x_1{x}_2{x}_3\right]}^T\hfill \end{array}$$

(3)

which can be used to formulate the model of the system as the following system of first order differential equations:

$$\left\{\begin{array}{ccc}{\dot{x}}_1\left(t\right)& =& x_2\left(t\right)\\ \\ {\dot{x}}_2\left(t\right)& =& -\frac{a}{2mb}exp\left(-\frac{x_1\left(t\right)}{b}\right)x_3^2\left(t\right)+g\\ \\ {\dot{x}}_3\left(t\right)& =& \frac{exp\left(\frac{x_1\left(t\right)}{p}\right)}{T_s}\left(k_{s}u\left(t\right)-i_s-x_3\left(t\right)\right)\end{array}\right.$$

(4)

All variables and parameters used in the model are described in

Table 1. Variables used in the mathematical description of a laboratory magnetic levitation system.

Variable Name	Descritpion	Value
$x_1\left(t\right)$	ball position	$x_1\left(0\right)=0.009\mathrm{m}$
$x_2\left(t\right)$	ball velocity	$x_2\left(0\right)=0$$\frac{\mathrm{m}}{\mathrm{s}}$
$x_3\left(t\right)$	coil current	$x_3\left(0\right)=0.1008\mathrm{A}$
$u\left(t\right)$	control voltage	$u\left(0\right)=0$
m	ball mass	$0.0571\mathrm{kg}$
g	gravitational acceleration	$9.81\frac{\mathrm{m}}{{\mathrm{s}}^2}$
$k_s$	gain of current driver	$2.5165\mathrm{A}$
$T_s$	time constant of the linear current driver	$0.031\mathrm{s}$
$i_s$	zero error of current driver	$0.1\mathrm{A}$
$a,b,p$	positive constants identified for the plant.	$a=0.017521\mathrm{H}$ $b=0.0058231\mathrm{m}$ $p=0.045626\mathrm{m}$

. We use exponential approximation for the electromagnetic strength formula. A detailed discussion of such choice can be found in [25].

It is a known fact that the linear controller can operate properly in the neighborhood of a chosen steady state. Performance of classical PID can be strongly improved, if the appropriate reference control value corresponding to a reference value is added to the generated control signal. The authors tested this solution with different types of non-integer PID controllers, see [26,27,28]. The constant value added to the control signal we can calculate by formula:

$$\begin{array}{ccc}\hfill u_r& =& \frac{-i_s+x_{3r}}{K_s}\hfill \end{array}$$

(5)

where

$$\begin{array}{c}\hfill x_{3r}=\sqrt{\frac{gmb}{a}exp\left(\frac{x_{1r}}{b}\right)}\end{array}$$

(6)

and $x_{1r}$ is the desired position of the ball.

The control loop of the described system was presented in Figure 2.

3. Method of Tuning of a Non-Integer Order PID Controller

One of the major tasks in the implementation’s field of the fractional PID controller is the formulation rules of tuning them. In this section, first we describe the concept of the PID controller with a fractional integrator, and then we propose two approaches of tuning it. The first one relies on introducing a fractional integral to a stabilizing PD controller and studying stability margins via the Nyquist criterion. The second one uses simulated annealing and multiple performance indicators. Both approaches are support methods requiring appropriate selection by the user.

3.1. Non-Integer PI${}^{\lambda }$D Controller

Podlubny (see [29]) proposed a generalization of the PID controller, namely the PI${}^{\lambda }$D${}^{\mu }$ controller, involving an integrator of order $\lambda$ and a differentiator of order $\mu$. As can be observed, when $\lambda =1$ and $\mu =1$, we obtain a classical PID controller, similar to when $\lambda =0$ and $\mu =1$ give PD, $\lambda =0$ and $\mu =0$ give P, $\lambda =1$ and $\mu =0$ give PI. All these classical types of PID are the particular cases of the fractional PI${}^{\lambda }$D${}^{\mu }$. However, the PI${}^{\lambda }$D${}^{\mu }$ is more flexible.

For laboratory active magnetic levitation systems, we used controllers by structured PI${}^{\lambda }$D. In the time domain, the equation for this controller’s output has the form:

$$\begin{array}{ccc}\hfill u\left(t\right)& =& K_{p}e\left(t\right)+K_i{}_0^C{\mathrm{D}}_t^{-\lambda }e\left(t\right)+K_d\dot{e}\left(t\right)\hfill \end{array}$$

(7)

where ${}_0^C{\mathrm{D}}_t^{-\lambda }$ is the Caputo derivative of order $-\lambda$ with respect to time, which is equivalent to the fractional integral of order $\lambda$, and:

• $K_p$ is proportional gain
• $K_i$ is integral gain
• $K_d$ is derivative gain
• $e\left(t\right)$ is control deviation in time t
• $\lambda >0$

In addition, the transfer function formula is given by the equation:

$$\begin{array}{ccc}\hfill G\left(s\right)& =& K_p+K_i{s}^{-\lambda }+K_{d}s\hfill \end{array}$$

(8)

3.2. Analytic Tuning Method with Nyquist Stability Criterion

This section describes the analytically procedure for determining the settings of the PI${}^{\lambda }$D${}^{\mu }$. They base the developed procedure on the Nyquist characteristic. The procedure of determining the PI${}^{\lambda }$D${}^{\mu }$ controller settings can be performed by the following steps:

• Determine the system work point.
• If the system is nonlinear, calculate the linearized system of the work point.
• Determine the stability region for the integer part of the controller. For example, if you want to create a controller PID${}^{\mu }$, you must calculate the stability region of a PI controller.
• Use a Nyquist criterion to determine the order of a non-integer order part of the regulator and this gain for a reasonable stability margin.

The advantage of this method is that we have an analytical description of the properties of the system—on the basis of which we can easily determine the stability area, gain, and phase margin. The most important disadvantage of this approach is that, for more complex systems, it can be difficult or not applicable.

3.3. Tuning by Simulated Annealing

Simulated annealing is a probabilistic technique for approximating the global optimum of a function. This algorithm uses a meta-heuristic to approximate global optimization in a large search space for an optimization problem—for problems where finding an approximate global optimum is more important than finding a precise local optimum in a fixed amount of time.

The advantage of this method is the ability to define any quality index for optimization, and there is no need to analytically determine the properties of the system. This method also allows the design of the regulator to bend known control characteristics. The disadvantage is the time needed to perform the optimization.

4. Design of the Controller for a Magnetic Levitation System

4.1. Construction PI${}^{\lambda }$D for Magnetic Levitation Systems Using the Analytic Method

To perform the calculations and simulations, the following assumptions have been made:

• $x_{1r}=10.8$ mm.
• Regulator must work for frequency bound between ${10}^{-4}$ and ${10}^3$ Hz.
• Order of approximation use for implemented non-integer elements must have value 5.

Based on assumptions, we can determine the work point of the system as:

$$\begin{array}{ccc}\hfill x_r& =& {[x_{1r},0,x_{3r}]}^{\mathrm{T}}\hfill \end{array}$$

(9)

In the next step, we calculate the transfer function of the linearized system in the work point:

$$\begin{array}{ccc}\hfill G\left(s\right)& =& \frac{vexp\left(x_{1r}/p\right)k_s}{T_s{s}^3+exp\left(x_{1r}/p\right)s^2-T_{s}ns-nexp\left(x_{1r}/p\right)}\hfill \end{array}$$

(10)

where $v,n$ have the form:

$$\begin{array}{ccc}\hfill v& =& -2\frac{a}{bm}exp\left(\frac{-x_{1r}}{b}\right)x_{3r}^2\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill n& =& \frac{a}{b^{2}m}exp\left(\frac{-x_{1r}}{b}\right)x_{3r}^2\hfill \end{array}$$

(12)

Then, we designed a classical PD controller which will stabilize the system. The stability region can by described by equation:

$$\begin{array}{ccc}\hfill K_p& <& \frac{n}{k_{s}p}\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill K_d& <& \frac{T_{s}n}{k_{s}p}exp\left(-x_{1r}/p\right)\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill K_d& <& T_s·K_{p}exp(-x_{1r}/p)\hfill \end{array}$$

(15)

Finally, a stability analysis based on the Nyquist criterion was performed. The Nyquist plot of the designed ideal (non-approximated) non-integer controller is presented in Figure 3 with a dashed line. Fractional controller parameters for these systems have value: $K_p=100$, $K_i=5$, $\lambda =0.75$ and $K_d=4$. As we can see, these parameters fully meet the assumptions.

4.2. Tuning of PI${}^{\lambda }$D Using Simulated Annealing

In this case, we can define the decision variables as: $K_p$, $K_i$, $K_d$, $\lambda$. The tests will be conducted for the following quality indicator:

• ISE: $\underset{0}{\overset{T}{\int }}{e}^2\left(t\right)\mathrm{d}t$
• ITSE: $\underset{0}{\overset{T}{\int }}t{e}^2\left(t\right)\mathrm{d}t$
• ISEDES: $\underset{0}{\overset{T}{\int }}\left(e^2\left(t\right)+{\dot{e}}^2\left(t\right)\right)\mathrm{d}t$
• IAE: $\underset{0}{\overset{T}{\int }}\left|e\left(t\right)\right|\mathrm{d}t$

where $e\left(t\right)=w_r-x_1\left(t\right)$. We use all of these qualities’ indicators because have different properties. ISE is one of the most used quality indicators. Since it is simple to use, it provides good results and in many cases can be calculated analytically. Its disadvantage is that it favors large errors, which means that, in the tuning process of the controller, errors at the beginning of the control process will have a greater impact on this process.

The ITSE quality indicator eliminates the problem of large fluctuations at the ends of the optimization horizon by multiplying the square of the control error by time. It is possible because the big value of errors at the beginning of the control process is leveled by multiplication by the small value of time; the same thing happens at the ending phase of the control process when the small value of errors are multiplied by the big value of time. This eliminates a significant influence of the error from the initial control phase on the value of the quality indicator.

ISEDES takes into account the rate of change of the error value, which allows for taking into account the dynamics in the optimization process. However, it has the same problems as quality indicator number 1.

IAE is presented, which treats all error values in the same way. However, its major disadvantage is that it can only be used in simulation environments because the analytical form of its solution does not exist. For more details, see [26,27,30,31].

Optimization start points have the following values:

• $K_p=500$
• $K_i=100$
• $K_d=6$
• $\lambda =0.5$

Simulation initial value:

• $x_1=9\mathrm{mm}$—since the system is inherently unstable, the starting point must be around the setpoint.
• $x_2=0\mathrm{m}/\mathrm{s}$—the sphere has no initial velocity.
• $x_3=0.1008$—value determined by Formula (6).

The optimal PI${}^{\lambda }$D settings for the system are collected in

Table 2. Result of tuning system using simulated annealing. The values of the $\lambda$ parameter are similar in all cases. The ITSE indicator sets the highest gain for proportional and integral parts of the regulator. In contrast, the IAE indicator sets much lower gains on these parts, but similar gains on the derivatives. It does not change the fact that the quality indicators assume values from a similar range. The ISE and ISEDES indicators set a similar gain of proportional and derivative parts. On the other hand, ISEDES is a weaker integrator gain and therefore gives the worst final value of the quality indicator.

Quality Indicator	${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	$\mathit{\lambda }$	${\mathit{K}}_{\mathit{d}}$	Quality Value
ISE: $\underset{0}{\overset{T}{\int }}{e}^2\left(t\right)\mathrm{d}t$	$493.168$	$109.581$	$0.998$	$0.822$	$4.967\times {10}^{-2}$
ITSE: $\underset{0}{\overset{T}{\int }}t{e}^2\left(t\right)\mathrm{d}t$	$523.684$	$181.68$	$0.959$	$2.606$	$2.116\times {10}^{-2}$
ISEDES: $\underset{0}{\overset{T}{\int }}\left(e^2\left(t\right)+{\dot{e}}^2\left(t\right)\right)\mathrm{d}t$	$499.257$	$74.654$	$0.982$	$0.875$	$7.209\times {10}^{-2}$
IAE: $\underset{0}{\overset{T}{\int }}\left|e\left(t\right)\right|\mathrm{d}t$	$455.075$	$146.7518$	$0.951$	$2.947$	$2.94\times {10}^{-3}$

. Position states of the magnetic levitation were shown in Figure 4. When analyzing the charts, we see that in all cases the system stabilizes around the set value in less than $0.1$ s. In all cases, you can see that there is a slight fix error. In the case of quality indicators 1 and 2, there is also a noticeable overshoot in the initial phase of the control process of magnetic levitation.

It is also visible in the values of the quality indicators collected in Table 2. Note that, for the quality indicators 1 and 3, the values of the $Kp$, $\lambda$, and $K_d$ parameters are similar. The highest values were also obtained for the indicators defined in this way. On the other hand, for the indicators defined with Formulas (2) and (4), the values of the $\lambda$ and $K_d$ parameters are similar. Additionally, it should be mentioned that the value of $\lambda$ is similar in all cases.

How can we see that the best result has been achieved when the quality indicator has form IAE? It should be noted that all quality indexes values are very close to each other. For this reason, the controller settings for the test in to laboratory system were selected values from IAE.

5. Realization Non-Integer Order PID Controller in a Real-Time Environment

5.1. Time Domain Oustaloup Approximation

The controller was implemented with the time domain Oustaloup method. This method was based on the classical Oustaloup approximation method, where fractional-order operator $G\left(s\right)=s^{\alpha }$ can be approximated by (see [32]):

$$\begin{array}{ccc}\hfill G_t\left(s\right)& =& K\prod _{i=1}^N\frac{s+{\omega }_i^{\prime }}{s+{\omega }_i}\hfill \end{array}$$

(16)

where:

$$\begin{array}{ccc}\hfill {\omega }_i^{\prime }& =& {\omega }_{min}{\omega }_u^{(2i-1-\alpha )/N}\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\omega }_i& =& {\omega }_{min}{\omega }_u^{(2i-1+\alpha )/N}\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill K& =& {\omega }_{max}^{\alpha }\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\omega }_u& =& \sqrt{\frac{{\omega }_{max}}{{\omega }_{min}}}\hfill \end{array}$$

(20)

The proposed approach is to realize every block of the transfer function (16) in the form of a state space system. Those first order systems will then be collected in a single matrix resulting in full matrix realization. This continuous system of differential equations will then be discretized. The time domain approximation has formulas:

$$\begin{array}{cc}\hfill A_k=& -{\omega }_k\hfill \\ \hfill B_k=& {\omega }_k^{\prime }-{\omega }_k\hfill \\ \hfill C_i=& 1\hfill \\ \hfill D_i=& 1\hfill \end{array}$$

This can be written in vector matrix notation

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}=& \left[\begin{array}{ccccc}{A}_1& 0& 0& \dots & 0\\ B_2& A_2& 0& \dots & 0\\ B_3& B_3& A_3& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ B_N& B_N& \dots & B_N& A_N\end{array}\right]\mathbf{x}+\left[\begin{array}{c}K{B}_1\\ K{B}_2\\ K{B}_3\\ ⋮\\ K{B}_N\end{array}\right]u\hfill \\ \hfill y=& \left[\begin{array}{ccccc}1& 1& \dots & 1& 1\end{array}\right]\mathbf{x}+Ku\hfill \end{array}$$

(21)

or in brief

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}=& \mathbf{Ax}+\mathbf{B}u\hfill \\ \hfill y=& \mathbf{Cx}+Du\hfill \end{array}$$

(22)

What can be immediately observed is that the matrix $\mathbf{A}$ is lower triangular. This is an extremely important in the case of this problem, as all its eigenvalues (poles of transfer function (16) are on its diagonal, so there is no need for eigenvalue products, which would lead to rounding errors. This is why discretization of (21) should have a structure preserving property.

The implementation of the algorithm requires the discretization of the control system designed in a continuous time domain. For preserving the stability attributes of the system in the discrete time domain, as supposed, the Tustin method also known as the bilinear transform for triangular matrix, has been chosen (see [33]). For more details about performance of time domain Oustaloup and realization fractional systems in real-time systems, see [28,34,35,36,37].

5.2. SoftFrac

The SoftFrac package was used to implement the regulators. SoftFrac gives Matlab user functionality to create a state-space system and transfer function fractional model of the non-integer dynamic element $s^{\gamma }$ based on approximations: Oustaloup (see [32]), time domain Oustaloup (see [38]), and LIRA (see [39]). The state-space fractional model class (ssf) inherits from standard state-space model class (ss) and have all of the functionality of this class, while the transfer function fractional model class (tsf) inherits from a transfer function model.

Because the SoftFrac library uses inheritance mechanism in implementation, the user can use all functionalities of parent classes ssf and tsf. In particular, the user can use this realization for Simulink simulation, system behavior analysis, and easy plotting of dynamic characteristics. In addition, we have added to the classes a construction method for easy conversion between those types, from ssf to tsf and vice versa.

6. Experiment Results

The controller performances for the magnetic levitation system has been investigated based on results from tuning methods in the following experiments:

• Introducing the ball into the magnetic field from below, near the working point by hand. Then, we waited for the ball’s position to stabilize at the working point.
• The experiment begins when the sphere is at the working point. It consists of knocking the sphere out of position by means of strokes.

The laboratory magnetic levitation system consists of: electromagnet, ferromagnetic sphere-position sensor, and current sensor power interface. The position measurement is done by the light source and sensor.

We can see in Figure 5 that the distance is measured from electromagnet to ball. The position measurement is based on the amount of light falling on the detector.

The composition of the laboratory stand on which the experiments were carried out included: PC, magnetic levitation plant (see [22]), and RT-DAC/PCI process board (see [40]) developed by INTECO. The computer on which the experiments were conducted used the Matlab 2015b version with the MATLAB/RT-CON real-time library installed. This computer used Windows 7 Professional, 64 GB RAM operating memory, and Intel Core I5 3.2 GHz.

RT-DAC/PCI I/O is a multifunction analog and digital I/O board dedicated to real-time data acquisition and control. The board contains a Xilinx FPGA chip. These boards can be reconfigured to introduce a new functionality of all inputs and outputs without any hardware modification. By default, the configuration version board accepts signals from systems and generates PWM outputs, and is equipped with the general purpose digital input/outputs (GPIO), A/D and D/A converters, timers, counters, frequency meters, and chronometers. The boards can be applied:

• using one of ready-to-use configurations distributed with the boards,
• using a new customer-specified configuration.

In this approach, the $n=5$ order of the non-integer order integrator approximation has been selected in the range $\omega \in \left[{10}^{-4},{10}^3\right]$ with sampling time of $0.001$ second.

Experiments have been conducted to validate the robust properties of the designed controller. The discrete Oustaloup approximation has been implemented in a real-time environment with the use of a RT-DAC process board and a MATLAB/RT-CON real-time library.

6.1. Analytic Method

In Figure 6a, one can find the stabilization position process in point 10mm. As can be seen, PI${}^{\lambda }$D controller stabilizes the system position in $1.5$ s.

The second experiment, presented in Figure 7a, shows the disturbance rejection for the controller. The system has been disrupted at about 2.5 s, 3.7 s, and 5.0 s.

6.2. Results of Tuning by Simulated Annealing

In Figure 6b, we see the stabilization position process for a regulator with settings from the tuning process. As can be seen, this controller needs more time to stabilize the system position, about 8 s.

The second experiment, presented in Figure 7a, shows the disturbance rejection for the controller. The system has been disrupted at about 2.5 s, 5 s, and 7 s. It should be noted that the regulators in both cases cancel the disturbance very well. This shows the high level of robustness of the regulators.

As shown in Figure 6, the regulator based on the SA method for the process of reaching the operating point has less overshoot than the regulator based on the analytical parameter selection criterion, while Figure 7 shows that both controls are designed to withstand disturbances in the same range.

7. Conclusions

Control design for the magnetic levitation system, while known, is still useful as a benchmark for both control quality and implementation purposes. The main issue that is known with fractional controllers is the problem of infinite memory. Direct implementation of fractional integrals (or derivatives) require either the summation of a series of signal samples with increasing length and varying coefficients or computing convolution integrals with singular kernels on expanding intervals. Both approaches are not practical, which is why approximations are required. In this paper, we have shown that implementing fractional controllers in real-time unstable systems is possible through fixed step realization based on time domain matrix representation. We have also studied the methods of tuning such controllers.

Initial attractiveness of optimization based approaches is unfortunately unsubstantiated by results. Minimizing the integral performance indicators can be realized by simulated annealing and, for example, Simulink created nonlinear models. The controllers obtained this way are, however, rather fragile, and problems with low frequency approximations of integrals lead to steady-state errors.

A more attractive approach is to use analytically based methods, such as the one we have proposed. In typical cases of such designs, the open loop system is stable, so one can use Bode plots for easy characterization. In the case of unstable systems like magnetic levitation, it is not possible. Here, we propose a different approach—designing a stabilizing integer order PD controller, and then introducing fractional integrals using a Nyquist criterion to keep stability margins. Such controller will of course have worse dynamics but should be able to compensate for incorrect identification and disturbance issues.

We are aware that our approach is not fully polished; however, it leads to a direct way for adopting fractional controllers by engineers—especially in the issue of robustness. While complicated equations of H${}^{\infty }$ designs might be difficult to adapt, frequency response shaping with fractional elements is much easier to understand.