Computation of Bifurcation Margins Based on Robust Control Concepts

. This article proposes a framework which allows the study of stability robustness of 4 equilibria of a nonlinear system in the face of parametric uncertainties from the point of view of 5 bifurcation theory. In this context, a branch of equilibria is stable if bifurcations (i.e. qualitative 6 changes of the steady-state solutions) do not occur as one or more bifurcation parameters are varied. 7 The work focuses speciﬁcally on Hopf bifurcations, where a stable branch of equilibria meets a 8 branch of periodic solutions. It is of practical interest to evaluate how the presence of uncertain 9 parameters in the system alters the result of analyses performed with respect to a nominal vector 10 ﬁeld. Note that in this article bifurcation parameters have a diﬀerent meaning than uncertain 11 parameters. To answer the question, the concept of robust bifurcation margins is proposed based 12 on the idea of describing the uncertain system in a Linear Fractional Transformation fashion. The 13 robust bifurcation margins can be interpreted as nonlinear analogs of the structural singular value, 14 or µ , which provides robust stability margins for linear time invariant systems. Their computation 15 is formulated as a nonlinear program aided by a continuation-based multi-start strategy to mitigate 16 the issue of local minima. Application of the framework is demonstrated on two case studies from 17 the power system and aerospace literature. 18

1. Introduction. Bifurcation analysis studies qualitative changes in the re-21 sponse of a nonlinear system (e.g. number and type of steady-state solutions) when 22 one or more parameters on which the dynamics depend are continuously varied 23 [27,21]. This is usually accomplished by selecting a few bifurcation parameters, typ-24 ically equal in number to the codimension of the studied bifurcation, based on their 25 importance for the system. This analysis approach is of recognized importance since 26 it allows complex dynamic behaviours to be characterized and an understanding of 27 the system to be gained. However, it does not provide indications on the robustness 28 of the results to uncertainties in the models. Let us consider for example the presence 29 of uncertain parameters allowed to vary within a prescribed range. These parameters 30 reflect the fact that uncertainty is ubiquitous in engineering systems and at any stage 31 of analysis (from preliminary to detailed). Unlike the bifurcation parameters, in prin-32 ciple they are not restricted in number (and are allowed to vary simultaneously) and 33 their influence on the dynamics may not be known a priori. It is then important to 34 estimate their effect, and in particular whether bifurcation points can move closer to 35 operating points deemed safe on the basis of analyses applied to the nominal system. 36 The study of robustness within a dynamical systems perspective can be attempted 37 by adopting singularity theory techniques (e.g. Lyapunov-Schmidt reduction) [18], as 38 shown by recently published research [20,8]. The central idea is to perform a reduc-39 tion of the original dynamics to a lower dimension map, whose singularities represent 40 transitions between qualitatively different bifurcation diagrams. Even though it is in 41 principle possible to track these singularities without computing explicitly the reduc-where x ∈ R nx and p ∈ R np are respectively the vectors of states and bifurcation 143 parameters, and f : R nx × R np → R nx is the vector field. In this work f is assumed 144 to gather smooth nonlinear functions (f ∈ C ∞ ). Therefore, the Jacobian matrix of 145 the vector field ∇ x f : R nx × R np → R nx×nx , denoted here by J, is always defined. 146 The vector x 0 is called a fixed point or equilibrium of (2. 1 J is singular at an equilibrium, i.e., it has a zero eigenvalue. The common feature 153 of static bifurcations is that branches of fixed points meet at the bifurcation point. 154 In the case of dynamic bifurcations, branches of fixed points and periodic solutions 155 meet. This case, also referred to as Hopf bifurcation, is the focus of this work and is 156 formally described by the following theorem.  (2.4) l 1 = 1 2ω H Re r, C(q, q,q) − 2B(q, A −1 B(q,q)) + B(q, (2iω H I n − A) −1 B(q, q)) .

186
The functions B : R nx × R nx → R nx and C : R nx × R nx × R nx → R nx are the tensors 187 of second and third order derivatives evaluated at x H , respectively. For example, for 188 vectors ξ, ς, χ ∈ R nx , B(ξ, ς) and C(ξ, ς, χ) are in R nx with components 189 (2.6) ∂x j x k x=x H ,p=p H ξ j ς k , i = 1, 2, ..., n x , ∂x j x k x l x=x H ,p=p H ξ j ς k χ l , i = 1, 2, ..., n x . at an initial point (x 0 , p 0 ), that there exist neighbourhoods X of x 0 and P of p 0 and 195 a function g : P → X such that f (x, p) = 0 has the unique solution x = g(p) in X.

196
Examples of numerical techniques to compute the implicit manifold g are  Raphson, arclength, and pseudo-arclength continuation [19], efficiently implemented 198 in freely available software, e.g., AUTO [14], and COCO [10]. 199 A general continuation problem, so called extended, can be formulated as follows and consider the restriction F (u, λ)| λ I =λ * I satisfying the IFT at some point (u * , λ * = 210 Ψ(u * )). Then F (u, λ)| λ I =λ * I defines a continuation problem for a d-manifold with If the model has uncertainties (which will be better characterized later), these can be 228 modelled with an operator ∆ u ∈ C q1×p1 with input z and output w. The effect of ∆ u 229 on M 22 can then be described by introducing the transfer matrices M 11 , M 12 and M 21 .

230
For example, in the case of parametric uncertainties, these will be simply static (gain) 231 matrices, while for the case of unmodelled dynamics these could also have dynamic 232 terms (e.g. low pass filters). The key point is that, by choosing these matrices, the 233 analyst can describe with a certain flexibility how the perturbation affects the nominal 234 system. Given this setting, Figure 1 shows the standard representation of LFT.

235
The central idea is thus to represent the uncertain system as a feedback of known 236 components (the transfer matrices M ij ) with uncertain (the operator ∆ u ) ones. In 237 practice, this is done by pulling out of the system the unknown parts, so that the 238 problem appears as a nominal system subject to an artificial feedback. Available 239 toolboxes [28] allow this operation to be efficiently performed and provide the analyst, 240 given a description of how the uncertainties affect the system, with the matrices M ij .

241
In order to formally define an LFT, let us denote by M ∈ C (p1+p2)×(q1+q2) the 242 partitioned transfer matrix (also termed coefficient matrix ) and let ∆ u ∈ C q1×p1 the uncertain operator. The LFT of M with respect to ∆ u is 245 defined as the map F : With reference to Fig. 1 where the uncertainties associated with n R real scalars δ i , n C complex scalars δ j , 261 and n D unstructured (or full) complex blocks ∆ D k are listed in diagonal format.

262
The identity matrices of dimension d i and d j take into account the fact that scalar 263 uncertainties might be repeated in ∆ u when the LFT  where s is the Laplace variable. Define now It can then be shown that F(M ν , ∆ ν ) = M 22 . This follows directly from where the diagonal structure of ∆ ν and the fact that 1 s = 0 have been exploited. (2.14)

301
A pictorial representation of the LFT F(M, ∆) defined by the operators in (2.14) is 302 given in Figure 2. The difference between the representations in Figure 1 and Figure 2, both describ-304 ing an uncertain system, is that in the former the system is described via its transfer 305 matrices, while in the latter a state-space representation is used. One can switch from 306 the first to the second representation by exploiting the fact that F(M ν , ∆ ν ) = M 22

308
The consequence of this change of representation is that the new block ∆ ν appears.

309
Correspondingly, the coefficient matrix M (2.14) now features the matrix M ν (2.12) 310 plus other matrices describing the effect of the uncertainties on the state-matrices.

311
Note indeed that the transfer matrices M 11 , M 12 and M 21 will also be expressed here 312 with their SS representation. Let us assume now that (2.11) is nominally stable (i.e.

313
A has all the eigenvalues in the left half-plane). Then the uncertain LTI system has 314 a purely imaginary eigenvalue if there exist ω > 0 (with s = iω) and a combination 315 of the uncertainties in ∆ u for which F(M, ∆) (2.14) is singular.

316
The advantage of this representation, which is key for the present work, is that The vector field depends now on δ, in addition to x and p. To highlight this, we 395 denote the uncertain vector field byf and the associated Jacobian byJ The objective of the work is then to compute the margins of stable equilibria from the 401 closest Hopf bifurcation for nonlinear systems affected by parametric uncertainties.

402
To better understand this, assume that the nominal system f has a Hopf bifurcation  It is often relevant to distinguish between supercritical and subcritical Hopf bifurca-419 tions, hence two distinct worst-case perturbations will be considered. For the sake 420 of readability, this distinction will be highlighted in the text when relevant but the 421 notation used will beδ in both cases.

422
In order to quantify the margin to the closest bifurcation, and thus to allow the 423 concept of worst-case uncertainty to be formalized, a metric for the magnitude of the perturbation must be adopted. The adopted metric should measure in some quanti-425 tative form the perturbation to which the system is subject. This task is arbitrary 426 and a common approach from robust control is followed [48] (see also section 2.2.2).
where (3.4a) exploits the property of interconnected LFTs, and ∆ u is a particular  The discussion above paves the way for the nonlinear program presented next, 494 which aims to compute the smallest perturbation for whichJ has a pair of purely 495 imaginary eigenvalues.
where X is the vector of optimization variables including: states x; uncertain param-499 eters δ; and frequency ω.X will indicate the solution vector gatheringx,δ, andω 500 respectively. Let us examine the constraints of the program. Eq. (3.5a) guarantees 501 that the solution (x,δ) corresponds to an equilibrium point for the system. Eq. (3.5b) 502 ensures thatJ has a pair of complex eigenvalues ν = ±ω, and Eq. (3.5c) bounds the 503 size of the perturbation matrix.

504
This is a similar optimization problem to that in (2.16), with two crucial dif-505 ferences: constraint (3.5a), and the addition of ∆ x in the block ∆ of F(MJ , ∆) (to 506 which, notably, constraint (3.5c) does not apply). Due to these differences, available 507 algorithms for µ cannot be applied to compute solutions of (3.5), thus alternative ways 508 should be pursued. Let us examine closely (3.5b), which prescribes singularity of an 509 LFT. According to the definition given in (2.9), necessary and sufficient condition As for (3.5c), this is a non-smooth constraint because of the maximum singular value 517 operator, but it can be drastically simplified by exploiting the structure of ∆ u (3.4c).

518
Indeed this constraint is equivalent to which is a set of linear inequalities in the optimization variables and the objective 521 function k m . Note that a similar relaxation would hold also for complex scalar uncer-522 tainties, not considered in this work.

523
Based on the previous discussion, the following smooth nonlinear optimization 524 problem is proposed to solve Program 3.1.
where n ctrs denotes the number of total constraints of the optimization.

528
The key idea behind Program 3.2 is to enforce singularity of the LFT (3.5b) by  This differs from the aforementioned works where the optimization was performed by 539 minimizing the nonsmooth functionσ(∆ u ). This is overcome here by considering 540 the relaxation commented in (3.6) and introducing the objective function k m as an 541 additional optimization variable.

542
Remark 3.1. Constraint (3.7b) consists of two (real and imaginary parts of the 543 determinant) nonlinear equality constraints in the variables X. By using Laplace 544 expansion of the determinant [1] and the fact that ∆ is structured, an analytical 545 expression for the gradient of (3.7b) with respect to δ and x can be obtained and 546 provided to the optimizer. As for ω, this is more tedious and therefore finite differences 547 are employed.

548
Note also that, from a continuation perspective, (3.7b) can be regarded as an analog  This is due to the fact that the frequency ω of the purely imaginary eigenvalues appear 554 explicitly in the constraint (and thus is an additional independent variable), which 555 is different from the test functions formulation. This is an important feature of the 556 developed approach, and possible ways to exploit it will be discussed later. can be written following (3.4) as where ν , x , and ω are unknown scalars described later. The following optimization 574 problem is then proposed to determine the worst-case perturbation for which both 575 conditions of the Hopf theorem are guaranteed to hold, that is, to calculate the margins 576 to the closest Hopf bifurcation point.

579
The first set of constraints (3.9a-3.9c) is identical to those in Program 3.2. The 580 constraints (3.9d-3.9e) instead ensure that the Jacobian linearized atp p has an ei-581 genvalue ν with real part ν (3.8d). Making use of a finite difference approximation, 582 it follows from the definition in (2.2) that l 0 = ν p . Therefore, existence of a solution 583 to Program 3.3 withˆ ν = 0 guarantees a Hopf bifurcation for the system. Hopf bifurcation to be computed.

639
Program 3.4. min X k m such that  Despite these measures, the risk of falling into local minima is still present.

705
In particular, the programs' initialization represents a critical aspect and thus a 706 continuation-based multi-start strategy is proposed. Assume that the optimizer has 707 found a solutionX to Program 3.2. The goal is then to provide the optimizer with a 708 set of initializations, derived fromX but possibly not leading the optimizer to find the 709 same solution, which allows an exhaustive optimization campaign to be performed.

710
The following extended continuation problem based on the constraints of Program 3.2 711 is first considered This can be recast in the formalism of (2.7) by setting 714 (3.13) introduced. The definition of g is arbitrary and various strategies can be pursued.

762
Crucially, the dimension is 1 irrespective of the number of uncertainties n δ , with the 763 drawback that these are now constrained to vary according to (3.14).    to these changes, two ordinary differential equations are added for V L and θ, and the 903 two algebraic equations become explicit equations for P L and Q L .

904
The resulting set of seven ordinary differential equations describing the power system,  θ)), where δ is the rotor angle, x T is the high side transformer reactance and x e is the 925 transmission line reactance.

926
The equations for the remaining three variables are not provided in [13]. The re-927 lationships for P L and Q L are derived here from the two aforementioned algebraic 928 equations in [7], which now allow an explicit expression for the load components to 929 be obtained since the phasor V L θ has a dedicated dynamic description (4.1f-4.1g).

930
As for E s , a relationship to the state variables is derived by considering the loadflow 931 equation for the circuit with the voltage source at the high side of the transformer.

932
This leads to:  Table 1 reports the values of the parameters used here for the power system  Table 1 Power system model parameters.    Table 2 for the parameters previously listed in Table 1. In order to show the connection between the sensitivity approach used in [13] and 982 the concept of robust bifurcation margin, a first type of analysis is discussed next. A 983 set of four parameters from the power system model is considered, namely x q , K A , 984 T A , and K f . Without loss of generality, only a subset of the parameters in Table   985 2 is selected to allow a more clear interpretation of the results. A subcritical value 986 of the loading level at which robustness of the plant is studied is then selected; this 987 is denotedl 0 according to the notation adopted in section 3.1. In all the analyses  Table 3 in terms of robust 1006 stability margin k m , frequencyω and worst-case perturbation for the normalized 1007 uncertainties.
1008 To further characterize this aspect, Figure 6 depicts the reciprocal of the robust 1063 bifurcation margin k m as a function of the frequency. The five curves represent the 1064 five cases considered in Table 4 and, unlike the one-shot tests discussed therein, are 1065 obtained by fixing the frequency in the optimization and computing the value of the 1066 margin at each frequency.

1141
The initial step to compute robust bifurcation margins is the definition of the   for this specific case, is known a priori to be the correct one. Moreover, at least for 1178 this case, Program 3.2 is able to detect the global minimum of the optimization.

1179
Another positive feature is that Program 3.2 has the frequency ω as a decision vari- that k m now achieves a smaller value than for c1. This is in accordance with the 1188 nominal analyses in Figure 7, for which c2 presented a smaller V H than c1. Thus, asV = 270 m s is closer to the nominal bifurcation speed for c2, a smaller robustness 1190 margin is expected. Note also that the two predicted frequenciesω are relatively close 1191 to the nominal ones. These interpretations thus give some confidence that an accurate 1192 prediction of the margin has also been obtained for c2.

1193
Other important information gathered in Table 5 is the type of closest Hopf bi-1194 furcations. Note that in order to obtain this result, the solver COCO was used to 1195 perform numerical continuation of the perturbed system, which also allowed verifi-1196 cation that the latter experienced a Hopf bifurcation atV 0 = 270 m s , as expected. only the case c2 will be considered.

1203
Uncertainties in two aerodynamic parameters are added to the set (4.9), namely, 1204 the terms A 012 and A 022 of the steady aerodynamics matrix A 0 (4.5). These corre-1205 spond to the lift and moment coefficients of the airfoil respectively, and are allowed 1206 to vary within ±20% from their nominal values. Table 6 shows the solutions provided  The supercritical case is consistent with the corresponding case in Table 5. Indeed,  Table 6. It is stressed that a quantitative interpreta-  Figure   1249 8 shows therefore that embedding the constraint on the Lyapunov coefficient in the 1250 bifurcation margin computation is successfully done by the optimization.

1251
The last part of the section is aimed at providing insights into the numerical as-1252 pects of the algorithms. As a preamble, it is observed that there are not definitive 1253 answers with respect to robustness to local minima or efficiency of the algorithms 1254 as these will depend on many aspects such as, for example, the type of vector field 1255 (not only size and degree, but also number of attractors) and the optimization al-1256 gorithm employed (which is an aspect that has not been investigated in this work).

1257
Investigation of these important features are left for future work.

1258
The execution time of Program 3.4 is larger than that of Program 3.2 (approxima-1259 tively 6s against 3s for the case with 7 uncertainties). Most importantly, the addition 1260 of the constraint on l 1 exacerbates the issue of local minima, especially when this is an 1261 active constraint. The set of strategies described in Section 3.3 were thus employed to 1262 obtain the results presented in Figure 8. Specifically, reinitializing the optimization