presents a singularity when S =W since in this case ϕ1 = ϕ3 while the trunk is always vertical (ϕ2 = 0). Therefore, here we rather concentrate on the control of the center of mass. There is a debate whether the center of mass position is directly controlled by the neural system or it is controlled indirectly through the body geometry [77]. Although this question has important neurological implications, from a mathematical point of view, indirect control of the center of mass for small displacements corresponds only to a reparametrization of the equations associated with direct control. The fact that many sensory organs for balance are in the head suggests that the control of the position of the head is also a possible control strategy, namely
C =Chead := S W
αhH
W −β− α W
, (3.8)
wherehH= 0.9(h−L)indicates the distance of the head from the middle of the hip.
As in [76], it is assumed that the delay and the control gains are invariant to system changes (see also [78]) and only the parameters of the mechanical model vary. The variations in the model parameters can be represented as uncertainties in the reduced inertiaI0, the gravitational termG0 and the configuration-dependent coefficient C. In order to concentrate the uncertain parameters, Eq. (3.1) is rearranged in the form
Ix(t) +¨ Gx(t) =−kpx(t−τ)−kdx(t˙ −τ), (3.9) where
I = I0
C and G= G0
C . (3.10)
Now all of the uncertainties are included in the normalized parameters I andG. The equation of motion can then be written in the first-order form
˙
x(t) =Ax(t) +BKx(t−τ) (3.11) with
x(t) = x(t)
˙ x(t)
, A=
0 1
−GI 0
, B= 0
1 I
, K=
−kp −kd
. (3.12)
We are interested in the robustness of this system against changes in the parametersI andG.
3.2 Robust stability analysis
A widely used tool for robust stability analysis is the pseudospectrum and the stability radius. The fundamental definitions of spectral value sets, complex and real stability radius were introduced in Chapter 1. In this section we redefine the problem for nonlinear matrix function (such as the quasi-polinomials of delay differential equations) and then extend the results of [76] further by considering structured perturbations following [12, 25, 29].
In [76], the robust stability of system (3.11) was determined using the unstructured complex stability radius associated with the perturbation of matrixA. If the system matrix is assumed to be uncertain, then the corresponding perturbed equation reads
x(t) = (A˙ + ˜A)x(t) +BKx(t−τ). (3.13)
24 CHAPTER 3 A human balancing model Substitution of the trial solutionx(t) =ceλt,c6=0into (3.13) yields
λI−(A+ ˜A)−BKe−λτ
c=0, (3.14)
which can be rearranged as
λI−A−BKe−λτ−1Ac˜ =c, (3.15)
ifdet λI−A−BKe−λτ
6
= 0. Taking the norm of both sides, using triangular inequality and simplification finally lead to
λI−A−BKe−λτ−1kA˜kkck ≥ kck
kA˜k ≥ λI−A−BKe−λτ−1−1
(3.16) Since (3.14) holds only if det λI−(A + ˜A)−BKe−λτ = 0, inequality in (3.16) gives the norm of the smallest perturbation matrix that satisfies this condition. This means that the perturbed characteristic function can have a zero on the imaginary axis if the norm of the perturbation matrix A˜ is larger than the right hand side of inequality (3.16). Since we seek the case when all of the eigenvalues remain on the left half complex plane, then the smallest perturbation matrix that satisfies this condition defines the complex stability radius. The formula is written as (see, e.g., [27])
rCA:=
sup
λ= iω ω≥0
λI−A−BKe−λτ−1
−1
. (3.17)
If the nominal system withA˜ =0is stable, then the perturbed system with anyA˜, kA˜k < rA
C is stable, too. Consequently, the robust stability boundaries for perturbations of different sizes are given by the contour curves ofrA
C.
In this unstructured case all of the elements of A are assumed to be perturbed by complex uncertainties (including the entries 0 and 1, which shall be invariant to system changes). This approach gives a conservative estimate of the actual robust stability associated with structured real perturbations. Furthermore, uncertainties of the mechanical parameters affect not only the system matrixA but also the input matrixB. This effect was also neglected in [76]. In this section, the robust stability analysis is presented in case of real-valued changes of the parametersIandG. Structured complex stability radius
The characteristic function of (3.11) can be written in the form D(λ) = det λI−A−BKe−λτ
=λ2I+G+kpe−λτ +kdλe−λτ. (3.18) In case of perturbations of the system parameters asI + ˜I andG+ ˜G, the weighted perturbation matrix can be introduced as [29]
∆:=
wII˜ wGG˜
, (3.19)
3.2 Robust stability analysis 25 wherewIandwGare the weights of the perturbations with respect to the nominal inertiaI and the nominal gravitational termG, respectively. If |I˜| < εI|I| and|G˜| < εG|G|, whereεI andεG are the radii of uncertainties, then the corresponding weights are
wI = 1
εII, wG= 1
εGG. (3.20)
For instance, if I is perturbed by maximum 2%, thenwI = (0.02I)−1. If wI → ∞orwG → ∞ then no perturbation on I or on G is allowed; cf [29]. This formalism allows the perturbations satisfying
k∆k2 =
I˜ εII
2
+
G˜ εGG
2
≤1. (3.21)
Thus, the allowed perturbations lie within an ellipse of main axesεIIandεGGin the plane( ˜I,G).˜ Following [29], the complex stability radius corresponding to the perturbation matrix∆can be calculated as
r∆
C =
sup
λ= iω ω≥0
kD(λ)−1w(λ)k
−1
, (3.22)
wherew(λ)is the complex weight function w(λ) =
λ2wI−1 wG−1
. (3.23)
If the nominal system is stable andrC∆ <1, then the system is robustly stable for any perturbations satisfying (3.21). Consequently, the robust stability boundary associated with the uncertainty radii εIandεGis given by the contour curver∆C = 1. Contour curvesrC∆ =zwith anyz ∈R+give the robust stability boundaries associated with the uncertainty radiizεI andzεG.
Structured real stability radius
Structured real stability radius can be defined using the real stability radius [19] and structured pseudospectrum [25], i.e.
r∆
R =
sup
λ=iωω≥0
υ D(λ)−1w(λ)
−1
, (3.24)
where υ(·) is defined by (1.29). Similarly to the previous cases, the robust stability boundaries associated with the uncertainty radii εI and εG is given by the contour curve r∆
R = 1, while the contour curvesr∆
R =z,z ∈R+ give the boundaries associated with the uncertainty radiizεI and zεG. Also note that in (3.24) the perturbations are limited to real parameters (such as the inertia) as opposed to (3.22), where complex perturbations are modeled. The complex stability radius is therefore significantly more conservative in general.
Robustness of balance control
Stability boundaries of the nominal system (without perturbation) in terms of the delayed feedback gains can be found in [57]. Robust stability boundaries in case of complex unstructured perturbation of the system matrixAwere provided in [76]. Here, robust stability analysis is presented for real-valued perturbation on the system parametersI andG. The stable regions of the nominal system
26 CHAPTER 3 A human balancing model in the parameter plane(kp, kd)were determined numerically using the semi-discretization method [1]. Stability radii and the robust stability boundaries are determined only in the stable regions.
In order to be able to compare the unstructured complex stability radius with the weighted structured real stability radius, we introduce the relative stability radius˜r. For complex perturbation, we use the ratio of the norm of the perturbation and the norm of the state matrix:
˜
rCA:= rAC
kAk. (3.25)
This number gives the relative complex perturbation that is allowed on the system matrixAwithout losing stability. For weighted structured real perturbation, we set the uncertainty radii toεI = 1 andεG= 1and introduce the relative stability radius as
˜
rR∆ :=r∆R εI=1
εG=1
. (3.26)
This number gives the relative real perturbation that is allowed on the inertiaIand the gravitational termGwithout losing stability. Now, the relative stability radiir˜A
C andr˜∆
R can directly be compared.
Figure 3.2 compares the robust stability boundaries and the pseudospectra calculated using complex unstructured (a-b) and real structured perturbations (c-d). Panel (a) shows the robust stability boundaries for complex unstructured perturbations of the system matrixAusing the same concept as in [76]. The equation under analysis is given in the first order form (3.11) and the robust stability boundaries are determined using the complex stability radius given by (3.17). Different contour curves associated with relative stability radiir˜A
C = 0,0.02, . . . ,0.08are presented. These contour curves correspond to0,2, . . . ,8% complex perturbation of the system matrixA.
Panel (b) in Fig. 3.2 illustrates the pseudospectrum associated with the control gainskp = 1140 and kd = 290. This point A is indicated by gray shading in panel (a). The three rightmost characteristic roots are indicated by small black dots and their pseudospectra associated with different perturbation levels are indicated by thin lines. The relative complex stability radius at this parameter point isr˜A
C = 0.041, i.e., the corresponding pseudospectrum just touches the imaginary axis. This indicates that there exists a complex perturbation matrix A˜ with kA˜k/kAk = 0.041 such that the perturbed system loses asymptotic stability. Thus, a4.1% complex perturbation ofA can already destabilize the system.
Panel (c) in Fig. 3.2 shows the robust stability boundaries obtained by the structured real stability radius according to (3.24) with εI = 1 and εG = 1. Different contour curves indicate different relative structured real stability radius r˜∆
R. For instance, the contour curves r˜∆
R = 0.1 indicates the robust stability boundaries corresponding to maximum 10% perturbation of the inertiaI and maximum10% perturbation of the gravitational termG, (see (3.21)).
Panel (d) in Fig. 3.2 shows the pseudospectra corresponding to the same control gains as in panel (b). This point B is indicated by gray shading in panel (c). Note, that A and B indicate the same control gains, but the pseudospectra are different. It can be seen that the change of the three rightmost characteristic roots for real perturbations is qualitatively different from that of the complex perturbations: the real characteristic exponents of the nominal system remain real and moves on the real axis as the perturbation changes. The wandering of the complex pair of eigenvalues is also different: they sharply drift to the imaginary axis in a specific direction. The relative real stability radius at this parameter point isr˜∆R = 0.44. This means that the system can lose stability if the real-valued perturbations satisfy
I˜ 0.44I
!2
+ G˜ 0.44G
!2
≥1. (3.27)
3.2 Robust stability analysis 27 Complex stability radius
Real stability radius
(a) (b)
(c) (d)
Stable
Stable
Unstable
Unstable
A
A
B
B
Pseudospectrum Pseudospectrum
Re(λ) Re(λ)
Im(λ)Im(λ)
kp
kp
0
0 0
0
0 0
0
0
500 500
1000 1000
1500 1500
2000 2000
kdkd 200
200 400
400 600
600 800
800
−15
−15
−10
−10
−10
−10
−5
−5
−5
−5
5
5 5
5
10 10
0.02 0.04 0.06 0.08
0.1 0.2 0.3 0.4 0.5
˜
rCA= 0 ˜rAC = 0.041
˜ r∆R = 0
˜
r∆R = 0.44
Fig. 3.2.Stability charts and corresponding pseudospectra with stance width ratioS/W = 1.2: (a) Complex stability radii; (b) Complex pseudospectra at parameter point A; (c) Real stability radii; (d) Real pseudospectra at parameter point B.
In other words, larger than 44% perturbation on the inertia I or on the gravitational termG can destabilize the system. This numerical example demonstrates that the complex stability radius gives a strongly conservative estimate of the relative stability radius. 4.1% relative complex perturbation on the system matrix A may already destabilize the system, while the size of the relative real perturbation onIandG, which gives an unstable system, is44%. The maximum relative stability radii in the robust stability diagrams show a similar tendency: r˜CA = 0.0814 whiler˜R∆ = 0.5801 (both points are indicated by dots in panels (a) and (c) of Fig. 3.2).
Note that the casesr˜A
C = 0andr˜∆
R = 0both give the conventional stability boundaries of (3.9) given by the D-curves
ifω = 0 : kp =−G, kd ∈R, (3.28)
ifω 6= 0 : kp = (ω2I−G) cos(ωτ), kd = ω2I−G
ω sin(ωτ). (3.29) The left bottom point of the stable region is given by
ωlim→0kp =−G (3.30)
ω→0limkd =−Gτ. (3.31)
The most robust point according to the real stability radius is located close to this point, i.e.,
kp &−Gandkd &−Gτ.
Figure 3.3 shows the effect of the application of different feedback signals. Panel (a) shows the case when the position and the velocity of the center of mass are applied as feedback signals according to [57, 76] (C = Ccom). Panel (b) presents the case when the feedback signal is the position and the velocity of the head (C =Chead). The uncertainty radii for both cases areεI = 0.2
28 CHAPTER 3 A human balancing model
(a) (b)
Unstable Unstable
Robust stable Robust stable
kp
kp
0 0 0
0 1000 2000 3000 1000 2000 3000
kd
kd
400 400
800 800
1200 1200
1.5 1.5
1 1
0.7 0.7
Fig. 3.3. Stable (white) and robustly stable (gray) domains for stance width ratiosS/W = 0.7,1,1.5for different feedback concepts: (a) position feedback of center of mass excursion; (b) position feedback of head excursion. The uncertainties of both the inertia and gravitational terms are20%.
andεG = 0.2, i.e., the maximum perturbations are20% for both parameters. Although the robust stability boundaries are different for the two feedback concepts, they are similar in topology, namely, there is no choice of(kp, kd), which can robustly stabilize the system for stance width ratiosS/W ranging from 0.7 to 1.5. This observation suggests that the nervous system might tune control gains to accommodate differentS/W configurations.
Figure 3.4 shows the change of the robust stability boundaries under different size of perturba-tions. The panels can be considered projections of the four-dimensional robust stability diagram in the parameter space(kp, kd, εI, εG). The uncertainty radiusεI for the inertia is kept constant in each row, while each column represents constant uncertainty radiusεGfor the gravitational term.
For instance, the caseεI = 20% andεG = 40% corresponds to perturbations which lie within the ellipse in the plane( ˜I,G)˜ defined by (4.23) with main axes0.2I and0.4G.
Stable regions for delayed feedback controllers are typically bounded by a straight vertical lines representing static loss of stability and curved boundaries representing dynamic loss of stability [6, 79]. It can be observed that the uncertainty of the inertia shifts the static stability boundary to the right, but does not affect significantly the dynamic stability boundary. The uncertainty in the gravitational term shifts drastically the dynamic stability boundary, but does not affect the static stability boundary. In case of the presented parameter combinations and stance widths, robust stability boundaries are more sensitive to uncertainties in the inertia than in the gravitational term.
Figure 3.5 shows the stabilizing feedback gains and the maximum relative stability radii for different stance width ratios. In panels (a,b) and (d,e), the stabilizing feedback gains are indicated by gray shading, the gains associated with the maximum relative stability radii are marked by solid lines. Panels (c) and (f) present the maximum relative stability radii as function of the stance width ratioS/W. Panels (a-c) shows the results obtained by unstructured complex stability radius according to [76], while panels (d-f) correspond to the structured real stability radius. For both types of perturbations the largest stability radii are associated with narrow stance widths.
There are two main differences between the complex and the real stability radii shown in panels (c) and (f). First, the maximum relative complex perturbation of the system matrix A without loosing stability is less than 10%, while this ratio for the real perturbation of the inertia and the gravitational term ranges between 48−63%. Thus, the allowed real perturbation on the actual mechanical parameters is much larger than the complex one. Second, the most robust derivative gainkdis located close to the lower stability boundary for the real perturbation, while it is about in the middle of the stable region for the complex perturbation. Note that the actual feedback gains fitted to human response to perturbation are in the lower left corner of the stable region [57], which