Graphical PID tuning method for uncertain fractional-order multivariable systems

1. Introduction

The industrial processes are often modeled by multivariate models. Fractional-order transfer function is always used to describe the industrial processes to improve the accuracy of multivariate system modeling. The parameters of models are usually uncertain due to the existence of measurement errors, machine aging and many other problems. Considering every input in the multivariable model may have impacts on outputs, the controller design of the multivariable model becomes very complex. Therefore, it is not feasible to treat the multivariable model as multiple univariate models directly. Meanwhile, the complexity of controller design for multivariable model will increase when the subsystems are fractional-order and the parameters of the multivariable system are uncertain. Multivariable system is confronted with decoupled issue when we design robust controller.

In past decades, many decoupling methods for the multivariable systems have been proposed. In literature [1], the SLC method is proposed to design the controllers of the multivariable system one by one in a certain matrix. Unfortunately, this method is very strict on the setting order of the controller. Li has proposed the practical multivariable control based on inverted decoupling and decentralized ADRC [2]. The modeling uncertainties, that significantly affect the robustness of the inverted decoupling control approach, are treated by ADRC outside of the decoupling structure. Jobrun has proposed a method for MIMO processes based on the multi-scale control scheme. The method tuning the controller parameters with the balanced sensitivity function [3, 4]. However, this method cannot be applied directly to the parameters uncertain multivariable system, because the uncertain parameters will bring huge calculation pressure to decoupling. There are also many optimization algorithms that can turn the parameter of controllers [5-7]. However, when the parameter of the system are uncertain, these algorithms will not be able to calculate the controller parameters. Therefore, a simplified decoupling method is needed to settle the coupling problem. The concept of an effective open-loop transfer function (EOTF) has been introduced and studied to design the controllers for multivariable systems by several researchers [8-13]. The core idea is to transform the multivariable system into a set of independent univariate systems. This method simplified the design of the controller while solving the coupling problem of the multivariable system. Vu [12] has proposed a method to design the controller and showed favorable control performance for multivariable system based on EOTF method and later studied by Jin [13].

Compared with integer-order calculus, the use of fractional calculus can be more accurate description of complex dynamic systems. For the fractional-order characteristic polynomials with time delay terms, the characteristic polynomial of fractional-order delay system is not a standard characteristic polynomial because of the existence of delay term and fractional order, while a quasi-characteristic polynomial that composed of the fractional power of the complex variable s and the delay term $e^{Ls}$. Therefore, the algebraic methods such as Routh-Hurwitz method, which is already used to test the stability of integer order linear control systems, cannot be directly applied to the stability criterion of fractional delay control systems. Of course, since the characteristic polynomial is a function of complex variables, the stability test principle of the existing integer order-delay system can be applied. The Kharitonov theorem has been proposed and studied to replace the parameter uncertain univariate system via multiple boundary equations for integer order linear control systems [14-17]. For each boundary equation, the region of the controller parameters satisfying certain robustness and stability can be determined. The intersection of all regions obtained by boundary equations is the controller parameter interval of the parameter uncertain univariate system. For fractional-order systems, Hamamci [18] proposed a method to find the stable region of fractional-order PID controller. Hamamci extended the D-factorization method for fractional-order PID controller to fractional-order domain and obtained the stable region of fractional PID. Using the existing stability theory of fractional calculus equation and Kharitonov boundary theory, Petráš et al. [19] proposed an algebraic method to test the stability of fractional linear system with uncertainties. Gao et al. [20] studied the relative stability of fractional-order systems and obtained the stable region of fractional-order PID controller. However, none of the above methods directly consider the uncertainty of the model parameters of the controlled object. Due to the working condition of the servo drive system, the model parameters also change with the system state. Therefore, the graphical tuning method should be based on the uncertain parameter model, otherwise, the set controller parameters will not meet the strong robustness requirements of the servo drive system, thus affecting the control accuracy. In response to this problem, Zheng et al. [21] proposed a graphical tuning method of fractional order proportional integral derivative (FOPID) controllers for fractional order uncertain system achieving robust D-stability. Hajiloo et al. [22] calculates the value range of the FOPID controller by multi-objective optimization. At the same time, Saidi et al. [23] gives the tuning method of the FOPID controller for fractional order uncertain system from the frequency domain. However, these methods are only for univariate systems.

In this paper, a new method is proposed to determine the robustness stable region of the parameter uncertain fractional-order multivariable system based on the EOTF method and Kharitonov method. In particular, the proposed FO-EOTF method is used to transform the parameter uncertain multivariable system into a set of independent parameter uncertain univariate systems. Then, the independent parameter uncertain univariate systems is numerical implementation by Oustaloup algorithm. Finally, the Kharitonov method is introduced to determine the robustness stable region of the parameter uncertain univariate system.

2. Effective open-loop fractional-order transfer function

The entire closed-loop of the system is shown in Fig. 1, where $R$, $U$ and $Y$ are the inputs, internal variables and outputs, respectively. $F_c$ is the PID controller. $F$ is the parameter uncertain fractional-order multivariable system. The model $F$ is shown as follow:

Fig. 1Closed-loop of the parameter uncertain fractional-order multivariable system

The parameters in model $F\left(s\right)$ are uncertain in a certain interval. The general expression is:

where $N_{ij}\left(s\right)$ and $D_{ij}\left(s\right)$ are the numerator terms and denominator terms of the transfer function $f_{ij}\left(s\right)$. ${\alpha }_{ij,m_{ij}}$ and ${\beta }_{ij,n_{ij}}$ are the fractional-orders of numerator terms and denominator terms. $a_{ij,h}$$\left(h=\mathrm{0,1},\cdots ,m_{ij}\right)$ and $b_{ij,k}$$\left(k=\mathrm{0,1},\cdots ,n_{ij}\right)$ are the coefficients of numerator terms and denominator terms, respectively, which satisfy ${\underset{\_}{a}}_{ij,h}\le a_{ij,h}\le {\stackrel{-}{a}}_{ij,h}$ and ${\underset{\_}{b}}_{ij,k}\le b_{ij,k}\le {\stackrel{-}{b}}_{ij,k}$. ${\tau }_{ij}$ is the time delays term of $f_{ij}\left(s\right)$ which satisfies ${\underset{\_}{\tau }}_{ij}\le {\tau }_{ij}\le {\stackrel{-}{\tau }}_{ij}$.

Fig. 2The concept of the FO-EOTF

The FO-EOTF method mainly utilizes the self-stabilizing ability of the system to solve the coupling problem, thereby reducing the difficulty of the coupling problem for the controller design. The FO-EOTF of loop $i$ is defined as the transfer function relating $u_i$ with $y_i$ where loop $i$ is open while all other loops are closed, as shown in Fig. 2, where $r_i$, $u_i$ and $y_i$ are the input, internal variable and output of loop $i$, respectively. ${\tilde{F}}_c^i$ is the remaining parts of $F_c$ without $f_{ci}$. ${\stackrel{-}{f}}^{ic}$ and ${\stackrel{-}{f}}^{ir}$ are the column $i$ and row $i$ of $F$ without $f_{ii}$, respectively. ${\stackrel{-}{F}}^i$ is the remaining parts of $F$ without column $i$ and row $i$. ${\stackrel{-}{r}}^i$, ${\stackrel{-}{u}}^i$ and ${\stackrel{-}{y}}^i$ are the remaining parts of $R$, $U$ and $Y$ without $r_i$, $u_i$ and $y_i$.

When ${\stackrel{-}{r}}^i=0$, ${\stackrel{-}{u}}^i$ can be written as:

Eq. (3) can be simplified as:

From Fig. 2, it is clear that $y_i$ is related to $u_i$ and ${\stackrel{-}{u}}^i$. The corresponding relationship can be written as:

When the frequency is below the cut-off frequency [10], the relationship between ${\stackrel{-}{F}}^i$ and ${\tilde{F}}_c^i$ is:

Eq. (5) can be written as:

Fig. 3The closed-loop of multivariable system after FO-EOTF

The influences from the controllers of the other loops are eliminated when tuning the controller of loop $i$. Each loop can be obtained by Eq. (7). Thus, the entire multivariable system $F\left(s\right)$ can be equivalent to a diagonal matrix system $F^{eff}\left(s\right)$, which is:

From Fig. 3, we can know that the structure of the whole closed-loop control loop after replacing the parameter uncertain multivariable system $F\left(s\right)$ with the parameter uncertain equivalent model $F^{eff}\left(s\right)$. The equivalent model can effectively reduce the interference caused by the coupling factors between the systems, and also greatly reduce the amount of computation. The equivalent transfer function $f_{ii}^{eff}\left(s\right)$ of loop $i$ in $F^{eff}\left(s\right)$ is only related to all elements in the system in Eq. (7). Assuming the general expression as:

where $h_{ii,*}\left(s\right)$ is consisted of $N_{ij}$ and $D_{ij}$ in Eq.(2), ${\tau }_{ii,*}$ is consisted of ${\tau }_{ij}$. Similarly, the FO-EOTF model similar to Eq. (9) can be obtained for each parameter uncertain multivariable system. However, the general expression of $f_{ii}^{eff}\left(s\right)$ in Eq. (9) is inconvenient to analyze and calculate because of the excessive time delay. Thus, Eq. (9) can be simplified as:

where $a_{ii,*}^{eff}$, $b_{ii,*}^{eff}$, ${\alpha }_{ii,m}$, ${\beta }_{ii,n}$ and ${\tau }_{ii}$ are all uncertain in a certain interval. In literature [9], it is pointed out that when the numerator and denominator coefficients of the transfer function are fixed, the larger the delay term is, the smaller the stability domain. First, considering the parameters of non-time delay term in $f_{ii}^{eff}\left(s\right)$ are fixed, the time delays are uncertain in a certain interval. Meanwhile, the coupling relationships among the coefficients ${\tau }_{ii,*}$ of the time delays are ignored. Then, the coefficients ${\tau }_{ii,a*}$ of the time delays in molecular term take the maximum value ${\stackrel{-}{\tau }}_{ii,a*}$. On the contrary, coefficients ${\tau }_{ii,b*}$ of the time delays in denominator term take the maximum value ${\underset{\_}{\tau }}_{ii,b*}$. Then, find the maximum values of ${\stackrel{-}{\tau }}_{ii,a*}$ and ${\underset{\_}{\tau }}_{ii,b*}$, that are $max\left\{{\stackrel{-}{\tau }}_{ii,a*}\right\}$ and $max\left\{{\underset{\_}{\tau }}_{ii,b*}\right\}$, and extract them out. Finally, the remaining parts after extraction are approximated by Taylor expansion. Take the 2×2 parameter uncertain system as an example:

where ${\underset{\_}{k}}_{*}\le k_{*}\le {\stackrel{-}{k}}_{*}$, ${\underset{\_}{T}}_{*}\le T_{*}\le {\stackrel{-}{T}}_{*}$ and ${\underset{\_}{\tau }}_{*}\le {\tau }_{*}\le {\stackrel{-}{\tau }}_{*}$. Through Eq. (7), the FO-EOTF of loop 1 can be written as:

where:

Under the premise of ignoring the coupling relation of the time delay, the fluctuation interval of the time delay can be obtained, that are ${\underset{\_}{\tau }}_{11}+{\underset{\_}{\tau }}_{22}\le {\tau }_{11,a1}\le {\stackrel{-}{\tau }}_{11}+{\stackrel{-}{\tau }}_{22}$, ${\underset{\_}{\tau }}_{12}+{\underset{\_}{\tau }}_{21}\le {\tau }_{11,a2}\le {\stackrel{-}{\tau }}_{12}+{\stackrel{-}{\tau }}_{21}$ and ${\underset{\_}{\tau }}_{22}\le {\tau }_{11,b1}\le {\stackrel{-}{\tau }}_{22}$.

Find out ${\stackrel{-}{\tau }}_{11,a1}$, ${\stackrel{-}{\tau }}_{11,a2}$ and ${\underset{\_}{\tau }}_{11,b1}$ according to the rules mentioned above. Then, $max\left\{{\stackrel{-}{\tau }}_{ii,a*}\right\}$ and $max\left\{{\underset{\_}{\tau }}_{ii,b*}\right\}$ can be obtained by comparison. For the convenience of display here, assuming ${\stackrel{-}{\tau }}_{12}+{\stackrel{-}{\tau }}_{21}\le {\stackrel{-}{\tau }}_{11}+{\stackrel{-}{\tau }}_{22}$, that means ${\stackrel{-}{\tau }}_{11,a2}\le {\stackrel{-}{\tau }}_{11,a1}$. Then, the remaining part of the extract is processed by Taylor expansion. Thus, Eq. (12) can be rewritten as:

In Eq. (13), the coefficient of the time delay is fixed, and the other coefficients are uncertain values in certain intervals. Due to the infinite dimensional properties of fractional-order operators, it cannot be directly numerically implemented. In this paper, the Oustaloup algorithm is used to approximate it by continuous transfer function. When using the Oustaloup algorithm to approximate fractional-order operator, the approximating frequency band and the order of the approximated transfer function need to select at first. Here, the approximating frequency band is assumed and the approximating order is $n$. The general transfer function of Oustaloup algorithm is expressed as:

where:

Therefore, Eq. (10) can be approximated by Oustaloup algorithm, which shown in Eq. (15) as:

3. Kharitonov theorem

It can be seen from Fig. 3, a multivariable system can be equivalent to multiple single variable systems after equivalent transformation. For general parameter uncertain univariate system, the corresponding boundary system can be obtained directly based on Kharitonov theorem. Then, the region of the controller parameters satisfying certain robustness for the parameter uncertain univariate system can be calculated. However, the premise of using the Kharitonov theorem is that the parameters of the uncertain model are linearly independent. The parameters $a_{ii,*}^{fo-eff}$ and $b_{ii,*}^{fo-eff}$ of the above FO-EOTF $f_{ii}^{fo-eff}\left(s\right)$ are determined by the coefficients in $f_{ij}\left(s\right)$, showing a linear correlation. In order to simplify the calculation, ignore the coupling between the parameters, readjust the intervals of the parameters as $a_{ii,*}^{fo-eff}\in \left[{\underset{\_}{a}}_{ii,*}^{fo-eff},{\stackrel{-}{a}}_{ii,*}^{fo-eff}\right]$ and $b_{ii,*}^{fo-eff}\in \left[{\underset{\_}{b}}_{ii,*}^{fo-eff},{\stackrel{-}{b}}_{ii,*}^{fo-eff}\right]$.

The controller used in this paper is PID controller. Taking loop $i$ as an example, the expression of the controller is:

where $K_{p,i}$, $K_{i,i}$ and $K_{d,i}$ are the proportions, integrals and differential terms of controller $f_{ci}\left(s\right)$, respectively. In order to satisfy the stability and certain robustness, a gain phase margin tester $A{e}^{-j\theta }$ is added to the loop $i$, as shown in Fig. 4. Thus, the characteristic polynomial of the loop $i$ can be obtained:

Fig. 4The closed-loop of multivariable system with gain phase margin tester

Substituting Eqs. (15) and (16) into Eq. (17), we get:

Eq. (18) can be written as:

Further algebraic manipulation of Eq. (19) results as:

where $q_{l,i}$ and $r_{l,i}$ are the real and imaginary parts of $s^l$ in ${\mathrm{\Delta }}_i\left(s\right)$ respectively, while $q_{l,i}$ and $r_{l,i}$ are uncertain but bounded. Thus, Eq. (19) can be rewritten as:

When Eq. (21) satisfies the Hurwitz stability, the entire uncertain system satisfies the stability condition. In order to further simplify the computation, the Kharitonov theorem is used to deal with the Eq. (21). The theorem indicates that Eq. (21) can be replaced by a set of boundary equations:

For each boundary equation ${\mathrm{\Delta }}_{\mathrm{*},i}$, there is a set of controllers $f_{c\mathrm{*},i}\left(s\right)$ that can make it stable. The intersection of controller parameters regions corresponding to each boundary equation is the final region of the controller parameters of the loop $i$, that is $f_{ci}\left(s\right)=\bigcap f_{c\mathrm{*},i}\left(s\right)$. Due to the existence of the gain phase margin tester $A{e}^{-j\theta }$, the controller $f_{ci}\left(s\right)$ can make the equivalent open-loop transfer function $f_{ii}^{fo-eff}\left(s\right)$ stable while satisfying certain robustness.

Let $s=j\omega$, the Eqs. (22)-(29) can be written as:

where ${\mathrm{\Delta }}_{R\mathrm{*},i}\left(\omega \right)$ and ${\mathrm{\Delta }}_{I\mathrm{*},i}\left(\omega \right)$ are the real and imaginary parts of ${\mathrm{\Delta }}_{\mathrm{*},i}\left(\omega \right)$, respectively. Let ${\mathrm{\Delta }}_{R\mathrm{*},i}\left(\omega \right)$ and ${\mathrm{\Delta }}_{I\mathrm{*},i}\left(\omega \right)$ equal to zero:

where $B_{R\mathrm{*},i}$, $B_{I\mathrm{*},i}$, $C_{R\mathrm{*},i}$, $C_{I\mathrm{*},i}$, $D_{R\mathrm{*},i}$ and $D_{I\mathrm{*},i}$ are the remaining parts after the common factor is extracted. Using the Cramer rule, the equations about $K_{p\mathrm{*},i}$ and $K_{i\mathrm{*},i}$ can be obtained:

When the parameters in the right part of Eqs. (33) and (34) are fixed, a boundary line can be calculated in the $K_{p,i}-K_{i,i}$ plane. The controller parameters corresponding to the point in this line can make the system in Fig.4 critical stable. Then, it is necessary to determine the regions of the controllers parameters by using the Jacobian matrix of Eqs. (33) and (34). The Jacobian matrix are:

When $J_{\mathrm{*},i}>0$, the left of the stable line facing the direction in which $\omega$ increases is the stable regions. When $J_{\mathrm{*},i}<0$, the right of the stable line facing the direction in which $\omega$ increases is the stable regions. In this way, the stability interval $f_{c\mathrm{*},i}\left(s\right)$ for each boundary equation can be determined. Thus, the final stable range $f_{ci}\left(s\right)=\bigcap f_{c\mathrm{*},i}\left(s\right)$ can be easily obtained.

4. Simulation

Assuming the model $F$ in Fig. 1 is a 2×2 parameter uncertain fractional-order multivariable:

where $K_{11}\in$[0.9, 1], $K_{12}\in$ [0.8, 0.9], $K_{21}\in$[1, 1.1], $K_{22}\in$[0.9, 1.1], $T_{11}\in$[1.7, 1.9], $T_{12}\in$[1.9, 2], $T_{21}\in$[1.5, 1.6], $T_{22}\in$[1.6, 1.7], ${\tau }_{11}\in$[0.45, 0.5], ${\tau }_{12}\in$[0.4, 0.45], ${\tau }_{21}\in$ [0.25, 0.3] and ${\tau }_{22}\in$[0.2, 0.25]. Assuming the loops need to satisfy the conditions of $GM\ge$ 8 dB and $PM\ge$ 30°. The final robust stable region of the loops can be calculated, as shown in Figs. 5-6.

Fig. 5Robust stable boundary lines of Δ1

Fig. 6Robust stable boundary lines of Δ2

Fig. 7Responses of controllers A1 and A2

Three sets of controllers are selected $A1$ ($K_p=$0.2, $K_i=$ 0.01, $K_d=$0.03), $B1\mathrm{}$ ($K_p=$0.2, $K_i=$ 0.08, $K_d=$0.04), $C1$ ($K_p=$1, $K_i=$0.5, $K_d=$0.1), $A2$ ($K_p=$0.25, $K_i=$0.02, $K_d=$0.03), $B2$ ($K_p=$ 0.25, $K_i=$0.05, $K_d=$0.06) and $C2$ ($K_p=$0.5, $K_i=$0.5, $K_d=$0.05) to verify the effectiveness of this method. Where, $A1$ and $A2$ are in robust stable region, $B1$ and $B2$ are in the stable region but outside of the robust stable region, $C1$ and $C2$ are outside the stable region. The responses of controllers$A1$ and $A2$ are shown in Fig. 7. The responses of controllers $B1$ and $B2$ are shown in Fig. 8. The responses of controllers $C1$ and $C2$ are shown in Fig. 9.

Fig. 8Responses of controllers B1 and B2

Fig. 9Responses of controllers C1 and C2

From Figs. 7-9, it clear that the controllers $A1$, $A2$, $B1$ and $B2$ can stable the system, but the controllers $C1$ and $C2$ cannot stable the system. In order to verify the effectiveness of the proposed method, the gain and phase margins of the controllers $A1$, $A2$, $B1$ and $B2$ are calculated to verify the robustness.

Fig. 10Bode plot of the controller A1 and A2

Fig. 11Bode plot of the controller B1 and B2

Fig. 10 shows the gain and phase margins of controller $A1$ and $A2$. The gain margin of controller $A1$ is in 28.7-31.6 dB, the phase margin is around 58 deg. The gain margin of controller $A2$ is in 32.9-36.7 dB, and the phase margin is around 51 deg. Fig. 11 shows the gain and phase margins of controller $B1$ and $B2$. The gain margin of controller $B1$ is around –18.4 dB, the phase margin is in 15.9-18.2 deg. The gain margin of controller $B2$ is around –21 dB, and the phase margin is in 32.5-37.1deg, which does not meet the limited conditions. Therefore, the controller in the robust stability region can make the system satisfy robustness.

5. Conclusions

In this paper, a unique graphical method are used to design the PID controller of the parameter uncertain fractional-order multivariable system. FO-EOTF is used to decouple multivariable systems and transform the multivariable system into a set of independent univariate systems. Based on Kharitonov theorem, a set of boundary systems is used to replace parametric uncertain fractional-order univariate systems. Then, the region of the PID controller parameters can be determined. In order to consider the robustness of the system, the gain phase margin tester is introduced. From the example, it can be seen that the proposed method can calculate the region of the PID controller parameters of the parameter uncertain fractional-order multivariable system. The parameters in the robust stable region can stabilize the system and meet certain robustness requirements. The parameters of the controller outside the robust stable region do not meet the requirements. The unique graphical approach helps to better understand how to determine the area of controller parameters.