Graphical Method for Robust Stability Analysis for Time Delay Systems: A Case of Study

Abstract

This chapter presents a tool for analysis of robust stability, consisting of a graphical method based on the construction of the value set of the characteristic equation of an interval plant that is obtained when the transfer function of the mathematical model is connected with a feedback controller. The main contribution presented here is the inclusion of the time delay in the mathematical model. The robust stability margin of the closed-loop system is calculated using the zero exclusion principle. This methodology converts the original analytic robust stability problem into a simplified problem consisting on a graphic examination; it is only necessary to observe if the value-set graph on the complex plane does not include the zero. A case of study of an internal combustion engine is treated, considering interval uncertainty and the time delay, which has been neglected in previous publications due to the increase in complexity of the analysis when this late is considered.

1. Introduction

The automatic control systems that are applied on the industrial process or on any production plant are previously analyzed, before they can be physically implemented. This analysis consists of obtaining some properties of interest such as qualitative and quantitative that are very helpful to achieve good performance in the closed-loop control system. The analysis is realized with the intention to predict the behavior that the control system will have when this is implemented. However, most of the time, this prediction lacks precision such that analysis is strongly based on the mathematical model of the physical process, and this does not represent in an exact way its dynamic behavior. This problem has been of great interest for scientific researchers in the last years, because it is desirable to get the best results when the control system is implemented but depends on these properties. A way of solving this problem is through the consideration of uncertainty on the mathematical model of the physical process. This uncertainty can be dynamic (see references [1–3]) or parametric (see references [4–6]). This chapter describes the qualitative property analysis of robust stability of a system. A case of study of an internal combustion engine is analyzed considering parametric uncertainty in the mathematical model. This case of study process is taken from references [5, 7], where conditions are obtained to verify the robust stability of the control system.

It is important to mention that in these research works several simplifications of the mathematical model were made, which may affect the real behavior of the system; one simplification performed by the authors is the cancellation of time delay in a part of the mathematical model, which has an influence that affects the stability property. The main contribution presented here is the consideration of the time delay on the mathematical model of an internal combustion engine, which is taken into account to obtain the property of robust stability. The methodology used is based on the application of the value-set concept to the particular case of the internal combustion engine (see reference [5]); to be more precise, it consists of the characterization of the value set of the resulting characteristic equation of the closed-loop system when the controller is connected. This controller assigns the poles in a position previously defined. Using this characteristic and applying the zero exclusion principle [8], it is possible to obtain robust stability conditions through a visual inspection of a graphic in the complex plane.

This chapter is organized as follows: Section 2 presents the elements used to verify the robust stability, Section 3 provides an abstract regarding obtaining the mathematical model of the internal combustion motor, in Section 4 the problem statement is set, in Section 5 the robust stability tools are used in simulation to verify robustness margin, and finally in Section 6 the conclusions of this work are presented.

2. Preliminary

In this section, we give the value-set concept, which is widely used to verify robust stability (see references [5, 9, 10]). The methodology that will be used consists of obtaining the value set for the characteristic equation, which results from the interval plants including time delay that are defined using the following definition:

Definition 1. A transfer function type interval plant is composed in the following manner:

$\text{g}\left(s,q,r\right)=\frac{n(s,q)}{d(s,r)}=\frac{{{\displaystyle \sum }}_{i=0}^m\left[q_i^{-},q_i^{+}\right]s^i}{s^n+{{\displaystyle \sum }}_{i=0}^{n-1}\left[r_i^{-},r_i^{+}\right]s^i}\forall q\in \mathcal{Q},r\in \mathcal{R}$E1

Now, if the interval plant includes the time delay, then it can be denoted by the following expression:

$\text{g}\left(s,q,r,e^{-\tau s}\right)=\text{g}\left(s,q,r\right)e^{-\tau s}\tau \in \left[0,{\tau }_{\max }\right]$E2

Notice that m < n, then g(s,q, r) is a set of strictly proper rational functions, Q and R are sets that represent the parametric uncertainty and are defined as follows:

$\begin{array}{l}\mathcal{R}\triangleq \left\{r={\left[r_1\cdots r_{n-1}\right]}^T:r_i^{-}\le r_i\le r_i^{+}\right\}\hfill \\ \mathcal{Q}\triangleq \left\{q={\left[q_1\cdots q_{n-1}\right]}^T:q_i^{-}\le q_i\le q_i^{+}\right\}\hfill \end{array}$E300

These kinds of sets are called as boxes. Equation (2) represents a type of systems known as interval plants. Both the numerator and the denominator of the transfer function have coefficients of uncertain values, which reside in a closed interval. The main interest in the study of this chapter is the robust stability analysis of feedback control systems, as the one depicted in Figure 1, whose main process is represented by an interval plant, with time delay and negative unitary gain feedback. The characteristic equation of this feedback system can be expressed in terms of the numerator, denominator, and time delay as follows:

$\text{p}(s,q,r,e^{-\tau s})=\text{d}(s,r)+\text{n}(s,q)e^{-\tau s}$E3

Figure 1.

Uncertain interval plant with time delay in negative feedback.

The term quasi-polynomial is used to describe these types of functions. Note that the above characteristic Eq. (3) is not a single equation but rather a family of an infinite number of characteristic equations. Then, this whole family must be considered if robust stability verification is carried out. This family is defined as follows:

$P_{\tau }\triangleq \left\{\text{p}\left(s,q,r,e^{-\tau s}\right)=q\in \mathcal{Q};r\in \mathcal{R};\tau \in \left[0,{\tau }_{\max }\right]\right\}$E4

The robust stability property is guaranteed if and only if the following equation is satisfied:

$\text{p}\left(s,q,r,e^{-\tau s}\right)\ne 0\forall s\in {\mathbb{C}}_{+}$E5

where ${\mathbb{C}}_{+}$ is used to denote the set of complex numbers with real part positive or equal to zero. From Eq. (3), it can be comprehended that the robust stability property of the dynamic system is very hard to verify using analytical methods, because this equation contains an endless number of equations. The goal of this work is to show a simple method to validate the robust stability property of time delay systems. The contribution of this work is based on the value-set characterization of the family of characteristic equations P_τ. The value set is defined here as follows:

Definition 2. The value set, denoted as V_τ (ω), of the characteristic equation of an interval polynomial with time delay is the set of complex numbers obtained by substituting s = jω in the polynomial:

$V_{\tau }\left(\omega \right)\triangleq \{p\left(s,q,r{e}^{-\tau s}\right):s=j\omega ,\omega \in \mathbb{R},q\in \mathcal{Q};r\in \mathcal{R};\tau \in \left[0,{\tau }_{\max }\right]\}$E6

where $\mathcal{Q}$ and $\mathcal{R}$ represent the set containing all possible values of the uncertainty of the parameters q_i and r_i expressed in a vector form as elements of the vector q and r. It is clear that the value set of P_τ is a set of complex numbers plotted on the complex plane for values of q_i, r_i, ω, and τ inside the defined boundaries. An important result that is applied in this chapter is the characterization of the value set for a characteristic equation as the one considered in the previous definition; this result is represented in references [11] and [12] with the following lemma:

Lemma 1. For each frequency ω, the value set V_τ(ω) is formed by octagons, where each one changes their geometry in function of the time delay τ. The coordinates in the complex plane of the vertices or corners of each octagon are given by the following formulas:

${\upsilon }_{i+1}=d_{i+1}\left(j\omega \right)+n_k\left(j\omega \right)e^{-j\omega \tau }$E7

${\upsilon }_{i+5}=d_{i+1}\left(j\omega \right)+n_h\left(j\omega \right)e^{-j\omega \tau }$E8

where

$i=0,1,2,3,$E1000

$k=\left(\gamma +i\right){\text{mod}}_4+1$E1100

$h=\left(\gamma +i+1\right){\text{mod}}_4+1$E1200

The term mod₄ represents the whole module base four operation, for example, mod₄ (3). γ can take integer values 0, 1, 2, and 3 depending on ωτ, see reference [9]:

$r\triangleq \{\begin{array}{cc}0& 2n\pi \le \omega \tau \le \frac{\pi }{2}+2n\pi \\ 1& \frac{\pi }{2}+2n\pi \le \omega \tau \le \pi +2n\pi \\ 2& \pi +2n\pi \le \omega \tau \le \frac{3\pi }{2}+2n\pi \\ 3& \frac{3\pi }{2}+2n\pi \le \omega \tau \le 2\pi +2n\pi \end{array}$E1300

The corresponding Kharitonov polynomials for the numerator and denominator of the interval plant are denoted by n_i (s) y d_i (s), respectively.

An example of a particular value set for fixed frequency and time delay is presented in Figure 2.

Figure 2.

Value set for a fixed frequency and time delay.

From the definition given for the value set, one can conclude that it includes all the values that the infinite family P_τ can take when it is evaluated in s = jω. Then, if the complex plane origin is contained in the value set V_τ(ω), this means that P_τ has roots located over the imaginary axis jω for some values of $\omega \in \mathbb{R}$. This, in fact, causes instability in the time delay feedback system. Consequently, the value-set technique can be used as an instrument that serves to validate the robust stability property. A question that arises in this point is how to show that P_τ does not have roots on the right half plane when a sweep over jω is done. The answer to this question is sustained in the following result known as the zero exclusion principle, see references [5, 11, 12]. This result will be applied to verify robust stability property.

Theorem 1: Suppose that p(s, q, r, e^‒τs) has at least one member stable for τ = 0, then p(s, q, r, e^‒τs) is robustly stable only if the value set satisfies:

$0\notin V_T\left(\omega \right)\forall \omega \ge 0$E9

Proof: See reference [8]

From the previous theorem, the robust stability problem is transformed into a problem where it only needs to verify that the value-set plot just avoids the zero of the complex plane.

3. Case of study mathematical model

In this section, the internal combustion motor is modeled. The idle mode operation condition is considered, which means the vehicle engine is running without accelerating; this is because the original objective is to reduce the fuel consumption when the vehicle is on and stopped, such as waiting for the green of the traffic light when circulating in the city. This fuel economy can be achieved increasing the air ratio of the fuel mixture, but this action causes instability in engine operation, resulting in a variation in the angular speed of the motor shaft. The mathematical model was compiled from article [7] and it is divided into three parts for better comprehension: (a) manifold chamber, (b) internal combustion chamber, and (c) rotational motion system.

3.1. Manifold chamber

The rate of chance of the pressure in the manifold chamber is affected by the current chamber pressure, the opening position of the throttle valve d (t), which controls the incoming air mass flow in a proportional way, and the outgoing flow that is proportional to the motor angular velocity n (t). The output of the chamber is the relative air pressure p (t). Then, the equation that gives the relationship between these two inputs and the output is given as the first-order differential equation as follows:

${\displaystyle \dot{p}}\left(t\right)+k_{2}p\left(t\right)=k_{1}d\left(t\right)-k_{3}n\left(t\right)$E10

where k₁, k₂, y k₃ are the respective proportionality constants. Applying the Laplace transform to previous equation, the following transfer function corresponding to the manifold chamber is obtained:

$p\left(s\right)=\frac{1}{s+k_2}[k_{1}d\left(s\right)-k_{3}n\left(s\right)]$E11

3.2. Internal combustion chamber

The combustion chamber produces the necessary torque to move the motor shaft. The torque generation subsystem can be modeled in a simple way having as inputs the spark advance (forward position of the rotor) a(t), the relative air pressure on the chamber p(t), the motor angular velocity n(t), and the fuel flow f (t). These variables contribute to the torque generation in a linear form. There is a time delay τ_d called induction-to-power-stroke delay, which affects the fuel control and the chamber pressure variables; this delay τ_d depends on the motor speed and the number of cylinders activated independently (denoted by n_c) as given by the following formula:

${\tau }_d=\frac{120}{n_{c}n(t)}$E12

Thus, the Laplace transform of the engine torque delivery T_e (s) generated for the combustion block is represented by the following equation:

$T_e\left(s\right)=e^{-{\tau }_{d}s}\left[k_{4}p\left(s\right)+k_{5}n\left(s\right)+k_{f}f\left(s\right)\right]+k_{6}a\left(s\right)$E13

Most of the time the delay is neglected, but when parameter uncertainty is considered, this has an influence on the control system stability, which is the reason this has been taken into account for analysis in this work.

3.3. Shaft rotational dynamics

Finally, the equations that represent the rotational can be obtained using Newton’s second law for rotational movement as follows:

$J{\displaystyle \dot{n}}\left(t\right)=T_e\left(t\right)-T_L\left(t\right)-k_{7}n\left(t\right)$E14

where J represents the rotational inertia, T_L is the torque of external load including its disturbances, and k₇ is an attenuation constant of the viscous friction that depends on the temperature, type of lubricant, and the gear that the motor uses. Therefore, the Laplace transform of the engine speed n (s) can be described as follows:

$n\left(s\right)=\frac{1}{Js+k_7}\left[T_e\left(t\right)-T_L\left(s\right)\right]$E15

Then, considering the transfer functions for each subsystem, the complete block diagram representing the mathematical model of the internal combustion motor can be obtained and implemented in Simulink MATLAB, for numerical simulation as shown in Figure 3.

Figure 3.

System block diagram in Simulink MATLAB.

The interest in this case of study is the relation between the variables d (s) and n (s) taking into account the time delay; then using block diagram algebra, the corresponding transfer function is given as follows:

$g\left(s\right)=\frac{k_1{k}_4{e}^{-{\tau }_{d}s}}{J{s}^2+\left(k_{2}J+k_7\right)s+k_2{k}_7+[k_3{k}_4-k_2{k}_5-k_{5}s]e^{-{\tau }_{d}s}}$E16

The previous transfer function will be considered as the mathematical model of the internal combustion motor.

4. Problem statement

It is important to mention that the engine mathematical model is expressed in terms of the parameters k_i and J. These parameters depend on the motor operating point; therefore, they cannot be considered as constant in the transfer function and they will be considered as uncertain terms that change depending on the operation point, but around the nominal parameters. The nominal values considered here are the same values as [7]: k₁ = 3.4329, k₂ = 0.1627, k₃ = 0.1139, k₄ =0.2539, k₅ = 1.7993, k₇ =1.8201, and the inertia J = 1. Making a change on the variables, the transfer function on the mathematical model can be simplified as follows:

$g\left(s\right)=\frac{a_1{e}^{-{\tau }_{d}s}}{s^2+a_{2}s+a_3+\left(a_{4}s+a_5\right)e^{-{\tau }_{d}s}}$E17

where the uncertain parameters are the ones included in a_i, the inertia J will be considered as a fixed number and equal to 1. Comparing the transfer functions (16) and (17), it can be observed that the new nominal parameters take the following values: a₁ = 0.8716, a₂ = 1.9828, a₃ = 0.2961, a₄ = 1.7993, a₅ = ‒0.2638.

A controller can be connected in the closed-loop system, as shown in Figure 4, with two objectives: the first one to regulate the angular speed and the second to improve the performance of the process. This controller is given by the following function:

$c\left(s\right)=\frac{50.0194s+26.3065}{s^2+9.8165s+33.1664}$E18

which assigns the system poles of the feedback control, considering a unitary feedback [13], at the following position: {−1, −2, −3, −4} on the complex plane. The controller was designed with the nominal parameters and time delay equal to zero. The characteristic equation of the closed-loop control system considering the uncertain parameters and the time delay has the following structure:

$p\left(s,q,r,e^{-T_{d}s}\right)=d\left(s,r\right)+n(s,q)e^{-{\tau }_{d}s}$E2400

where

$d\left(s,r\right)=s^4+\left(a_2+9.82\right)s^3+a_3+9.82a_2+33.17)s^2+\left(9.82a_3+33.17a_2\right)s+33.17a_3$E19

$n\left(s,q\right)=-a_4{s}^3+\left(a_5-9.82a_4\right)s^2+(9.82a_5-33.17a_4+50.02a_1)s+33.17a_5+26.31a_1$E20

The problem considered in this work is to determine the property of robust stability of the internal combustion motor when uncertainty is included in the new nominal parameters a_i. This property is directly related to the characteristic equation.

Figure 4.

Control for angular velocity of the internal combustion motor.

5. Robust stability verification

First, the characteristic equation is considered without uncertainty, taking into account the nominal parameters and the time delay, which means:

$p\left(s,q,e^{-T_{d}s}\right)=s^4+r_1{s}^3+r_2{s}^2+r_{3}s+r_4+\left(q_1{s}^3+q_2{s}^2+q_{3}s+q_4\right)e^{-{\tau }_{d}s}$E2700

where r₁ = 11.80, r₂ = 52.93, r₃ = 68.67, r₄ = 9.82, q₁ = −1.80, q₂ = −17.93, q₃ = −18.67, q₄ = 14.18.

Using Lemma 1, the corresponding value set for $\omega \in [0,.5]\text{y}{\tau }_d\in \left[0,9.9\right]$ is obtained, which is represented on the graphic shown in Figure 5.

Figure 5.

Value set without uncertainty in the parameters.

From the previous graphic, it can be appreciated that the value set does not reach the zero of the complex plain, but is very close; thus, it can be determined that the maximum delay the control system supports to preserve the stability is 9.9 s. Several numerical simulations can be run for different time delay limits to appreciate more clearly if the value set reaches the origin, for instance, considering the time delays ${\tau }_1=2.5,{\tau }_2=9.9,{\tau }_3=13$, which correspond to the stable, oscillatory, and unstable system responses, respectively, as shown in Figures 6–8, respectively.

Figure 6.

Transient response for stable behavior.

Figure 7.

Transient response for oscillatory behavior.

Figure 8.

Transient response for unstable behavior.

Finally, considering an uncertainty of 10% in each of the parameters, the following characteristic equation is obtained, where an overestimation was made to obtain interval polynomials with time delay:

$p\left(s,q,r,e^{-T_{d}s}\right)=s^4+r_1{s}^3+r_2{s}^2+r_{3}s+r_4+\left(q_1{s}^3+q_2{s}^2+q_{3}s+q_4\right)e^{-{\tau }_{d}s}$E2800

where r₁ = [11.60,11.99],, r₂ = [50.95, 54.90], r₃ = [61.80,75.53], r₄ = [8.84,10.80], q₁ = −1.98, −1.62], q₂ = [−19.72, −16.13], q₃ = −29.25, −8.08], q₄ = [11.01,17.35].

Once again, using the lemma, the value set of the characteristic equation is obtained now with uncertainty in the coefficients; This can be appreciated in the graphic of the figure 9.

Figure 9.

Value set with uncertainty in the parameters.

The previous value set was obtained for a value range of $\omega \in [0,1.5]$, ${\tau }_d\in \left[0,2.8\right]$ and like the previous case it can be observed that it barely reaches the zero of the complex plain; therefore, it can be assumed that the maximum time delay the control system supports is 2.8 s. Therefore, it can be clearly appreciated that parameter uncertainty decreases the delay margin that the control system supports and it is important to take this into account when the robust stability property is being verified. It is necessary to note that due to the overestimation about the parameters to represent the characteristic equation as a polynomial delay interval, the results obtained only warranty enough robust stability conditions; nevertheless, the next simulation shows how instability is presented for each of the values contained in the uncertainty, on the control system, when having a delay of τ_d = 4.1 s (see Figure 10).

Figure 10.

Step response with uncertainty in both parameters and time delay.

6. Discussion and conclusion of results

A robust stability verification methodology based on a graphic method was mentioned in this chapter. First, the mathematical model with nominal approximated parameters must be obtained; then, a control algorithm connected in feedback with the plant is analyzed when in the plant of the closed-loop system, uncertainty is introduced in both parameters and time delay. The main tools used are the value-set graph and the zero exclusion principle, then the analytical problem to verify the robust stability is transformed into a graphical problem, which consist of checking if the value-set plot for frequency interval does not include the origin. The case of study of an internal combustion motor was presented. The time delay appears when modeling the corresponding section of the engine combustion chamber of the motor. It is important to note that by considering the uncertainty in the parameters of the mathematical model, the time delay influences more in the control stability system; so, it can become significant even when it is very small and therefore the inclusion in the analysis is important. Then, with the tools presented here, a robustness margin for parameters and time delay can be obtained.