5.3 Numerical Stability Analysis
5.3.1 Numerical Calculation of Stability Diagrams
In [HajduInsperger16] stability diagrams are analytically constructed for continuous sys-tems with delayed state feedback and dead time compensation using a Smith Predictor.
The analytical constructions are based on the D-subdivision method and Stepan’s formulas, as presented in [St´ep´an89]. The D-subdivision method allows to calculate the stability boundaries by separating real and imaginary part of the characteristic equation. Moreover, Stepan’s formula determines the number of unstable characteristic exponents in each domain of the stability chart, see [St´ep´an89].
In the following, this method of illustrating stability regions is transferred to the analysis of digital control systems. Therefore, the analytical construction of the stability charts are replaced by a numerically calculation of eigenvalues of the discrete system augmented by delayed states as derived in Eq. (5.34) and Eq. (5.37). Similar numerical methods, as semi-discretization, are presented in [InspergerSt´ep´an02] and [InspergerSt´ep´an11]. The analytical and numerical stability charts are calculated and compared in the following for a mass-spring system.
Example: Mass-Spring System With Delayed Feedback
The objective is to calculate numerically a stability chart as considered in [HajduInsperger16], where the chart is derived analytically. The system matrix of the mass-spring system is given as
A=
"
0 1
−a 0
#
, (5.39)
with a = 0.5 and hence, the system matrix is stable. The input matrix is given as B =
"
0
−1
#
. (5.40)
A state feedback is used and defined as
u=kTx(t−τ), (5.41)
with delayed state x(t−τ), time delay τ and control gain k= [p, d]T ∈R2. Hence, the closed loop system dynamics reads
˙
x=Ax(t) +BkTx(t−τ). (5.42) The analytically calculated stability chart from [HajduInsperger16] for the continuous system (5.42) is shown on the left in Fig. 5.6 for time delay τ = 1 and varying control gain p∈[−1,1] and d∈[−1,2]. The stable regions are shaded in gray and the system is only stable for proportional gain p ≥ −a. In the following, the steps for the numerical calculation of the stability charts are given.
1. Discretization
The continuous system is discretized using Eq. (5.13) and Eq. (5.14). The sampling time T0 = 0.05s is chosen. Then, the discrete system matrix and input vector read
Ad=
"
0.99 0.05
−0.025 0.99
#
, Bd=
"
−0.001
−0.05
#
. (5.43)
-1 -a 0 1
-1 0 1 2
p
d
ω
ω= 0
-1 -a 0 1
−1 0 1 2
p
d
Figure 5.6: Comparison of analytical stability chart according to [HajduInsperger16] (left) and numerically calculated chart (right).
5.3 Numerical Stability Analysis 105
2. Augmentation Time Delay
The discrete system is augmented by delayed states. The number of delayed steps is given as ntotal = τtotalT
0 = 20, see Eq. (5.26). Hence, the augmented system order yields na = 2 + 2·20 = 42, as given in (5.28). The augmented state vector at time point k reads xa,k = [x1,k, x2,k, x1,k−1, x2,k−1, . . . , x1,k−20, x2,k−20]T ∈R42. (5.44) The augmented discrete matrices are
Ad,a=
"
Ad 02,40 I40,40 040,2
#
, Bd,a=
"
Bd 040,1
#
, (5.45)
as derived in (5.31) and (5.32).
3. Delayed Augmented Feedback
The state feedback controller is stated with delayed discrete states. Using the augmented delayed state, the delayed discrete state feedback controller reads
uk =kTxk−20. (5.46)
Then, the control gain vector kT is augmented to
kTa = [01,next, p, d], (5.47) with zero matrix 01,next and dimension 1×next. It follows that the augmented feedback controller can be rewritten to a function of the augmented state vector from Eq. (5.44) with
uk=kTaxa,k. (5.48)
Finally, the delayed discrete closed loop system reads
xa,k+1 =Ad,axa,k+Bd,akTaxa,k = Ad,a+Bd,akaT
| {z }
=Ad,a,cl
xa,k, (5.49)
with discrete augmented closed loop system matrix Ad,a,cl.
4. Evaluation Eigenvalues
Now, the eigenvalues of the closed loop system matrix Ad,a,cl can be evaluated over a grid of values of control gains p, d. The closed loop system is unstable, if there exists a discrete eigenvalue zk,l of Ad,a,cl for a parameter constellation pk, dl, with
|zk,l|>1. (5.50)
Otherwise the closed loop system is stable and the point pk, dl is shaded grey in the stability chart to mark a stable region.
The stability chart on the right in Fig. 5.6 shows the numerical evaluation over a grid of pk∈[−1,1] with step size 0.01 and dl ∈[−1,2] with step size 0.025. Hence, the discrete augmented closed loop system matrix Ad,a,cl is evaluated 24321 times.
The figure demonstrates that analytical (left) and numerical (right) stability charts are similar and therefore the numerical approach can generate stability diagrams with sufficient accuracy. However, the numerical approach cannot provide exact D-curves, which represent stability boundaries in the analytical chart. But on the other hand, the numerical method has other advantages. The stability charts can be easy calculated using numerically methods and further investigations on the eigenvalues are possible. In the following investigations regarding the maximum magnitude and damping ratio, are derived.
The maximum magnitude |zmax| of the discrete eigenvalues zi with i ∈ {1,2, . . . , na} represents how far away the eigenvalues of the stable systems are from stability boundary.
It is calculated by
|zmax|= max{|z1|,|z2|, . . . ,|zna|}. (5.51) Furthermore, the dominant damping ratio ξdom of the most dominant eigenvalue of the stable system is calculated, since it determines in large part the damping behavior. The most dominant eigenvalue is given as the continuous eigenvalue with the greatest stable real part. Therefore, the transformation from discrete to continuous eigenvalues is necessary.
The transformation is given as
λi = ln(zi) T0
, (5.52)
as for instance discussed in [Lunze16].
Then, the dominant damping ratio ξdom can be calculated using Eq. (3.38).
Figure 5.7 shows the stability charts of the closed loop delayed system with maximum absolute eigenvalue on the left side and dominant damping ratio on the right side. The red regions are favored since they mark low maximum magnitudes of the discrete eigenvalues or high damping, respectively. Both favored properties can be considered separately. The regions with the largest intersection between stability and damping is desired. For instance a control gain with p=−0.25, d= 0.4 is a good choice for stability and damping.
Furthermore, time simulations using the continuous time system with delayed feedback from Eq. (5.42) are performed with various control gains in order to verify the calculated stability charts. The system is deflected by the initial condition x0 = [1,0]T. Figure 5.8 illustrates time simulations with the following control gains, chosen from the stability charts:
5.3 Numerical Stability Analysis 107
Figure 5.7: Numerically calculated stability charts with color coded maximum absolute eigenvalue of the discrete eigenvalues (left) and color coded damping ratio of the dominant eigenvalue (right).
• a)p=−0.25, d= 0.4, small maximum magnitude and high damping
• b)p=−0.15, d= 0.4, small maximum magnitude, average damping
• c)p= 0.25, d= 0.4, maximum magnitude 1, no damping, marginally stable
As the stability charts show, the simulation with control gain a) has the best damping behavior such that the steady-state x1,ss = 0 can be reached fast. The simulation using control gain b) has average damping behavior and the simulation using gain c) has no damping and the system is stable but not asymptotically stable, as calculated in the stability chart.
Figure 5.8: Time simulation of the first state x1 of the continuous delayed closed loop mass-spring system with various control gains to check the stability charts.