- Email: [email protected]
Definition and Computation of a Nonlinearity Measure
Copyright © IFAC Nonlinear Control Systems Design, Tahoe City, California. USA, 1995
Definition and Computation of a Nonlinearity Measure Frank Allgower* *Inlttitut Jir SYlttemdyno.milc und Rege/unglttechnilc, Univerltitiit Stuttgo.rt, 70550 Stuttgo.rt, FRG, phone: +49-711 - 6856193; Jo.x: +49-711-6856371; emo.il: o.lIgower~rult . uni-lttuttgo.rt.de
Abstract. Nonlinearity measures are a means to quantify the size of nonlinearity in the 1/0behavior of nonlinear systems. This paper is concerned with the definition and the computation of nonlinearity measures. We review some previous results on the computation of nonlinearity measures and also demonstrate a novel approach to efficiently compute nonlinearity measures based on convex optimization. An example illustrates the theory and demonstrates the complexity of the computations involved. Key Words. nonlinearity measures; quantification of nonlinearities; convex optimization
1. Introduction
discrete time systems that can be represented by second order Volterra series. In (Nikolaou, 1994) a different measure is proposed that quantifies the average nonlinearity of a nonlinear I/O-system in discrete time. A completely different nonlinearity measure related to the achievable closed loop performance is described in (Allgower, 1994b). A comprehensive discussion of different measures can be found in (Allgower , 1994a). The practicality of the different definitions depends of course on the application and goal at hand , but also on the computational complexity of its calculation . In this paper the focus is on the derivation of applicable schemes for the computation of a particular nonlinearity measure that appears to be of practical usefulness. We will review our previous results on the calculation of this nonlinearity measure and introduce a new method for its computation based on convex optimization. The paper is structured as follows: In Section 2 the nonlinearity measure is introduced and discussed. Section 3 is concerned with the computation of this nonlinearity measure . In Section 3.1 calculation of a lower bound based on the computation of Fourier coefficients is described . Section 3.2 describes a novel approach based on convex optimization. A short discussion about possible applications of nonlinearity measures in connection with nonlinear control system analysis and design together with an example is given in Section 4.
Not all nonlinear systems have the same "degree" of nonlinearity. There are systems with nonlinearities that have no significant influence on the I/O-behavior. Other nonlinear systems have a behavior that is very different from the behavior of any linear system. By definition , both types of systems are nonlinear. However the first type of systems could be termed "mildly" nonlinear, whereas the second type could be called "strongly" nonlinear. For analysis and synthesis purposes it is often crucial to distinguish between them. In many cases a distinction simply based on intuition is not sufficient. The introduction of nonlinearity measures is an attempt to systematically approach this problem by quantifying the size of nonlinearity in the I/O-behavior of nonlinear systems. The term nonlinearity measure was first used in (Desoer and Wang , 1981) . Without the need for an explicit computation, it was shown there that certain linear feedback laws have a linearizing effect on nonlinear systems. The first systematic application of nonlinearity measures to the analysis of nonlinear control systems and the synthesis of nonlinear controllers was reported in (Allgower and Gilles, 1992) together with algorithms for their computation . Independently nonlinearity measures were also applied in (Nikolaou, 1994) and (Sourlas and Manousiouthakis, 1992) where also computational schemes for their calculation are developed. The measure considered in (Sourlas and Manousiouthakis , 1992) is equivalent to the measure defined in (Desoer and Wang, 1981) and to the one used in this paper. The computational scheme presented in (Sourlas and Manousiouthakis, 1992) is however restricted to
2. Definition of a nonlinearity measure We consider general nonlinear (i.e. not necessarily linear) , multivariable I/O-systems with input u and output y . The input signals u are
257
assumed to belong to the normed space U a of admissible input signals. The output signals y are elements of the normed space of output signals y. The I/O-behavior is described by the nonlinear operator N that maps input signals u E Ua into output signals y E Y :
Other sensible definitions for non linearity measures can be defined as well. For example a normalization of (3) by the "size" of y (= N[u)) may be useful in some cases (Allgower, 1994b): ~u N
y = N[u).
(1)
lim yet)
=0
for
u(t)
==
0
"It.
Definition 1 The nonlinearity measure
sup
In general the computation of nonlinearity measure
(2)
GEy UEU
IIC[u]- N[u]llpy , lIull puG
(3)
3.1. A lower bound on
The "size" of the nonlinearity of an I/O-system is thus defined as the normalized largest difference between the output of nonlinear system Nand a linear system C measured by the norm 11 . Ilpy in y. The difference is taken with respect to the worst-case input in U and with respect to the linear system C* that approximates N best. An example for a specific choice of nonlinearity measure is given by restricting U to consist of L 2 -signals with bounded energy lIullL 2 :S u and the norms 11 . IIpy and 11 . IIpu G both being the L 2-norm . The set U contains the subset of inputs that we consider relevant in the controlled loop. It is easy to see that the nonlinearity measure
Restriction ofU to a suitable set can greatly facilitate the computation of the nonlinearity measure. For reasons of clearer presentation we restrict ourselves in this subsection to SISO systems with continuous inputs and outputs u, y E C[O, 00) satisfying
with C : Ua ........ Y being a linear I/O-operator, 1I·lIpy a norm in y , 11·llpua a norm in U a , and U being the set of inputs considered.
IIC*[u]- N[u]lIpy = 0
(5)
3. Computation of the nonlinearity measure
of a BIBO-stable, causal I/O-system N : U a ........ Y for inputs signals u E U ~ Ua is defined by the nonnegative number ",U . f 'l'N= 10
IIG[u) - N[u)llpy IIN[u]lIpy .
In the remainder of this paper we restrict our attention to the nonlinearity measure defined by (3).
We consider static as well as dynamic operators N, but restrict our attention to causal, BIBOstable systems in this paper. This setup describes a very general class of systems that comprise finite dimensional systems in continuous and discrete time as well as infinite dimensional systems. Without loss of generality, we assume t-+oo
. f
= J~y ;~ft
limT-+oo lIullvT limT-+oo lIyllvT
< <
00
(6)
00
with
lIyllVT =
liT
-
T
y2(t) dt.
(7)
0
We consider the special class of sine-inputs u(t) = A . sin(wt),
(8)
where the amplitude A and frequency ware taken from the sets
A EA = {A E R+ 10:S A :S Amar} (9) wEn = {wER+ IWmin:SW:SWmar}. This specific set of inputs, that we will denote by U$ in the sequel, is thus parametrized by A and w :
(4)
U.
= {u 1 u(t) = A ·sin(wt),A E A,w En} .
(10)
For linear systems the dynamic behavior is completely described by the response to sine-inputs with one arbitrary amplitude and w E R+. Despite the additional variation of the amplitude, the behavior of a nonlinear system N is of course not completely described by the response to input signals in U•. Any nonlinearity measure based on
The behavior of nonlinear system N for the input signals considered is thus the same as that of the linear system C* and therefore N behaves linear in a certain sense. We will refer to this as: "N is linear in U". It must be stressed that the value of the nonlinearity measure depends of course on the choice of input signals U considered .
258
U. instead of U will only be a lower bound on the measure we are really interested in. In particular we suggest the following approximate nonlinearity measure X~' for N:
r
IIG[u] - N[u]lIvT lIullvT
If we consider the set of input signals U. instead of the single input (8), then by the parametrization of U. the relation to be proved follows immediately. 0
(15)
Calculation of the lower bound (11) thus only involves computation of Fourier coefficients for different periodic output signals. For simple nonlinearities, like for example static memory less nonlinearities, an analytical calculation is possible in analogy to the methods for determining describing functions (e.g. Gelb and Vander Velde (1968)). For systems that are more complex, the Fourier coefficients, and thus the lower bound, can be calculated in an efficient and reliable way by numerical techniques (e.g. Natke (1988)). The numerical computation of X~' is only based on the simulated output signal. Thus no restrictions have to be assumed on the type of nonlinear systems that can be considered. Therefore for any system that can be simulated, including systems described by differential-algebraic equations, systems with nonsmooth nonlinearities , etc. , this lower bound on the nonlinearity measure can be calculated. From (11) it follows immediately that if N is a linear I/O-operator, the value of the lower bound is also equal to zero. From X~' = 0 it can however only be concluded that for each input signal u E U there exists a linear system Gu (depending on the respective u), so that the outputs of Nand Gu do not differ for this specific input signal u.
Theorem 1 If the set U. is given by (1 0) and the measure X~' is defined as in eq. (11), then
3.2. Numerical computation of the nonlinearity measure through convex optimization
U, A
XN =
.
f
:~E, J~gT:'rr;.,
(11)
Note that sup and inf are interchanged in (11) as compared to (3). This, together with U. ~ U, results in X~' being a lower bound on 4>#
< A-.U _'I'N,
U,
XN
where
1I ·lIp»
and
1I ·llpu
a
in
(12)
4># are also defined as (13)
The advantage of this particular lower bound is that it can be calculated in a very efficient way. Under some technical assumptions, the response of nonlinear system N to a sine-input (8) with fixed A E A and fixed wEn is the sum of two functions (Isidori and Byrnes, 1989): (i) the periodic steady-state response YIO , and (ii) the transient response Ytr that decays to zero, because of the assumed stability of N : YN
= N[u] =
Ytr
+ YIO ·
(14)
The periodic signal YIO can thus be expanded in a Fourier series with amplitude coefficients Ak 00
YIO
= Ao +
2: A
k'
sin(kwt
+ 4>10) .
10=1
U
XN'
= sup AEA wEn
1 -.
A
2( A) Ao w,
~ AHw , A)
+~
2
In this section we want to demonstrate an alternative method to approximately compute the nonlinearity measure 4>}j. defined by (3) . This approach is based on solving an appropriate nonlinear programming problem numerically. In order to apply numerical methods we have to make two approximations:
' (16)
10=2
where Ak(W , A) are the amplitude coefficients of the ph harmonics of the steady-state response YlO to inputs u(t) = A . sin(wt) .
sketch of proof: As a first step we consider only one input signal ofform (8) with fixed amplitude A and fixed frequency w . From the principle of harmonic balance (e.g. Vidyasagar (1978)) we know that the output YL of the linear system G* that approximates N best in the sense of eq.(l1) is given by the first harmonic of the output of N yL(t)
= G*[u](t) = AI ' sin(wt + 4>1)'
• We restrict the set of inputs considered to a finite set, that we will call Uc in the sequel. • We predefine a fixed structure for the unknown best linear approximation that only contains a finite number of variable real parameters. We will make no further restriction on Uc neither with respect to the type of inputs considered (sinefunctions, steps, random inputs, etc.) nor with respect to the number of input signals. The problem we want to solve is: Find G* E g, such that
(17)
with AI , 4>1 as in (15). Using this G* , it follows after some straightforward algebraic manipulations that inf lim IIG[u] - N[u]llvT GEgT-oo
(18)
= lim IIG*[u] - N[u]lIvT = T-oo
< 259
where 9 is the space of linear operators mapping u E Ua into y .
very satisfying results. Further theoretical investigations on the appropriate choice of basis functions and with respect to convergence properties are however needed. Calculation of the approximate solution B~c of the nonlinearity measure ,p via numerical optimization as described has some desired properties due to the convexity of the problem (Nemirovsky and Yudin, 1983):
Theorem 2 The infinite dimensional optimization problem (19) is a convex programming problem. proof: We have to show that
f :9
~
#-
R+ with
f(G) = max IIG[u]- N[u]lIpy UEUc lI u ll pU4
(20)
is a convex function with respect to G and that 9 is a convex set . As 9 is a Banach space, convexity of 9 follows immediately. In order to show the convexity offunction f we define another function fft with ii E Ue fixed and
fft(G)
= IIG[ii] -
N[ii]llpy·
(21)
This function is convex with respect to G by the definition of convexity: 11
A Gdii] + (1-A)G 2[ii]- N[ii]lIpy (22) = IIAGdii] + (1-A)G2[ii]- AN[ii] - (l-A)N[ii]lIpy ::; AIIG1[ii] - N[ii]lIpy + (1-A)IIG 2 [ii] - N[ii]lIpy-
Because the maximum of a finite number of convex functions (scaled by constants lIullpU4) is again a convex function (Ekeland and Temam, 1976), f(G) = maxftEUc lIulfpU4 . fft(G) is also a convex function. 0
Therefore 8~c can be calculated numerically in an effective and reliable way. If the set of inputs Ue considered for calculation of B~c is not dense in U , then B~c will only be a lower bound on ,p#-. Further investigations will concern the use of an appropriate choice of basis functions for Ue , such that U e is dense in U. The computation of 8~c is more involved than computation of the lower bound X~·. If only sineinputs u E Ue = U. (compare (10» are considered and the norm in (11) and (19) is equally defined for both cases , then X~· ::; B~· . Thus B~' is a tighter bound on ,p~' than X~· . Another advantage of the proposed numerical solution via convex optimization as compared to calculation of X~· is that we are not restricted to consider only sineinputs. Any input signals that appear to be appropriate for the problem at hand can be used for the computation of B~c.
The convex infinite dimensional optimization problem (19) can be approximated through a convex finite dimensional problem by a suitable convex parametrization of the space of linear operators. A possible parametrization for discrete time SISO systems is for example given by assuming a moving average structure (step response model) for linear system G : n
G[u](k) =
LXi ·u(k - i) ,
(23)
;=0
where the scalars Xi are the step response coefficients to be determined . This choice preserves the convexity between the variable vector of parameters Z = [Xl , .. . , xn]T and function fin (20). For continuous time SISO systems for instance a convex parametrization
4. Example
n
G[u] =
LXi ·G;[u],
• There is a unique global optimum. No further local optima exist. • The cost to determine this global optimum increases the most polynomially with the number of variable parameters. • The cost grows the most linearly with respect to the number of input signals u contained in Ue . • In order to determine the global optimum, the response of nonlinear system N to each input u E U e has to be computed only once. Thus a calculation for complex systems is as feasible as a calculation for a simple system. • There exist very powerful algorithms to solve convex optimization problems (e.g. (Nesterov and Nemirovsky, 1993; Grotschel et al. , 1988» .
(24)
Possible applications of nonlinearity measures for the analysis of nonlinear control systems are of course manifold and include control structure selection, validation of nonlinear models and calculation of the best linear approximation. Besides their importance for a general analysis of nonlinear systems, nonlinearity measures are especially interesting in connection with approximate 1/0linearization where they can be used as a tool for a quantitative assessment of the size of nonlinear-
;=1
with G;[u] fixed stable linear systems , is also possible. The extension to multivariable systems is straightforward. It is clear that the quality of the approximate solution depends on the choice of approximation functions. Our own experience and experience collected in a different context (Boyd et al. , 1988; Scherer , 1994) shows however that such finite dimensional parametrizations lead to
260
To demonstrate the flexibility of the approach, we base the computation of eJ:je on a discretized model of (25) that is derived using finite forward differences with relatively fast sampling 6T = 0.2. Considering a set of inputs of 24 different stepfunctions with amplitudes between -0.6 and +0.6
ity that remains after approximate linearization. In this context, nonlinearity measures can also be used to determine the number of approximation terms needed in order to achieve a satisfying linearity. For a detailed discussion see (Allgower, 1995; Allgower, 19946). In this section we want to quickly compare the computational methods presented using a simple example. This will also serve to demonstrate the numerical effort needed and to discuss implications of the choice of input signals U . We consider the following nonlinear SISO system with unstable zero dynamics Y=
Uc1
for the computation of eJ:je leads to a value of
(30) for the uncompensated system. h(t) in (29) denotes the unit step function . The convex optimization was carried out using the MINIMAX function of the MATLAB Optimization Toolbox (Grace, 1992). For this system, with Uc1 as in (29) , a simulation horizon of 100 timesteps, and a step response model with n = 25 coefficients, about ten seconds CPU time were needed on the same VAX-Station as above. It is also possible to compute e in continuous time. For this we consider 24 pulse functions of different amplitude
(25)
with scalar input u, scalar output y, and state :I: = [Xl , X2JT . We are interested in the behavior of
:- ········:··.:.:11' ::.::.:.·:.::.::.:..:::.: .. ::::·:1[·.·:':".::"'::.::','··············· .....·1j I
:
i
i
...................
I
,
I
,
i.................,,.
et
i
.. _-_ ... ,;
!
··· .... ·: .. ·.. ········· .. ·t· .. ·· .. ··· .. ·· .. ·, ··· ···· .. ·...... \
!
with
A E {-0.6, -0.55, .. . , +0.55, +0.6} } (29)
Xl
Xl = -3XI + 4X2 + 2xr + u X2 = -Xl + X2 - 0.5x~,
= {ulu(t) = A . h(t)
u = 0.2
.. ...... !
u
~--+--+--- i
u
= 0.4
= 0.6
Uc 2
= {ulu(t) = A . h(t) -
A . h(t - 50) with
A E {-0.6, -0 .55 , ... , +0.55, +0.6}} Fig. 1. Step responses of uncompensated systems (25) .
and use the parametrization of the best linear approximation according to (24) with only five fixed, stable linear systems
system (25) in a neighborhood around the steady state :1:$ = [O,OV and assume that the input is constrained by
lul:::; 0.6
C;[u]
(26)
.
xt·
eJ:je2 = 0.67.
(33)
Computation time was about four seconds. Using convex optimization it is also possible to calculate 1 with respect to the input signals described by (27). The resulting measure has a value of
et·
= {u I u(t) = A· sin(wt) with A E {0.3, 0.6} and w E {10, l,O.5,O.l,O .05 , 0.01} }
X·
r. ' ,T; E {.5,.99,1.73,2.82,4.45} . (32) s ; +1
Optimization was carried out on the same VAX computer as above using ACSL and an (almost) standard SQP algorithm, leading to a value for the nonlinearity measure of
System (25) displays strongly nonlinear 1/0behavior as demonstrated by the step-responses shown in Figure 1. 1 Calculation of the lower bound of the nonlinearity measure for N using a set of twelve sine-inputs U. I with different amplitudes and frequencIes U. I
(31)
(27)
eJ:j·1 = 0.62
leads to a value of
(34)
for the continuous as well as for the discrete setup . The resulting best linear approximations C* of N derived by all three methods, virtually show identicall/O-behavior. From these figures it can be seen that the methods to compute the nonlinearity measure presented in Section 3.1 and Section 3.2 lead to very similar values of the nonlinearity measure. The 1 value of the lower bound only differs from eJ:jel and eJ:je2 by less than ten percent, sustaining our experience that the lower bound is usually rather tight. It can also be seen from the above that the nonlinearity measure is rather insensitive
(28) for system (25) . The value of the nonlinearity measure, i.e. the gain from input u to the difference between the output of nonlinear system N and the output of its best linear approximation, is thus almost as big as the gain of the nonlinear system itself. This confirms our observation that system N is strongly nonlinear. Computation of 1 for system (25) using the set of inputs (27) in the simulation environment ACSL (Mit, 1991) only needed fractions of a second on a VAX-Station 4090.
xt·
xt·
261
to the choice of input signals as long as the relevant operating regime and dynamic range is captured properly. Even the type of parametrization for the best linear approximation doesn't seem to alter the results in a significant way. System (25) can be approximately 1/0linearized by feedback u = -1.9xr
+ 0.05x~ -
0.14xIX2
+v
(35)
despite its unstable zero dynamics (Allgower, 1995). With this feedback the non linearity measure of the compensated system is
OJj·l =
X~;l = 0.05 .
(36)
Figure 2 shows step responses of the compensated system (25) with feedback (35). The considerably "more linear" behaviour of the compensated system is clearly reflected by the value of the nonlinearity measure. .......... : ........ .
....... :
. ..... .1
..... 1. ......... ..i ...............L · i !
j
= =
v -0 .6 v = -0.4 v -0 .2
v = 0.2 v = 0.4 v = 0.6
Fig. 2. Step responses of system (25) with compensation (35) .
5. Conclusions In this paper we discussed a measure for quantification of the size of nonlinearity in the 1/0behavior of a large class of dynamical systems. This nonlinearity measure captures the difference in the I/O-behavior between the nonlinear system and its best linear approximation. Computation of this nonlinearity measure is usually rather involved. We have given an overview over two methods to approximately calculate the nonlinearity measure: a lower bound based on evaluation of Fourier coefficients, and a novel approach using convex optimization . These approaches are computationally very attractive. Even for rather complex systems of high order a numerical solution is possible. Acknowledgement: Fruitful discussions with C. Scherer are greatfully acknowledged. We also want to thank D. Odenthal, W . Schobbert, and M. Klander for their invaluable help in implementing the computational codes.
6. REFERENCES Allgower, F . (1994a). Entwurfvon naherungsweise E/ A-linearisierenden Reglern fur nichtlineare Systeme. Interner Bericht. ISR, Universitat Stuttgart. Allgower, F. (1994b). NichtlinearitatsmaBe - Ein Werkzeug zur Analyse und Synthese ni(:htlinearer Regelkreise. In : Entwurf Nichtlinearer
262
Regelungen (S. Engell, Ed.). Oldenbourg Verlag. Allgower, F. (1995). Definition and computation of a non linearity measure and application to approximate I/O-linearization. Internal Report 95-1. ISR, Universitat Stuttgart. Allgower, F. and E.D. Gilles (1992). Approximate input/output-linearization of nonlinear systems. AIChE Annual Meeting, Miami, FL. Boyd, S., V. Balakrishnan, C. Barratt , N. Khraishi, X. Li, D. Meyer and S. Norman (1988). A new CAD method and associated architectures for linear controllers. IEEE Trans. Automat. Contr. AC-33(3), 268-283. Desoer, C.A. and Y.-T. Wang (1981). Foundations of feedback theory for nonlinear dynamical systems. IEEE Trans. Circ. Syst. CAS27(2), 104-123. Ekeland, I. and R. Temam (1976). Convex analysis and variational problems. North-Holland . Amsterdam . Gelb, A. and W . E. Vander Velde (1968) . Multiinput Describing Functions and Nonlinear System Design. McGraw-Hill. New York, NY. Grace, A. (1992). Optimization Toolbox for use with MATLAB. The MathWorks, Inc .. South Natick, MA. Grotschel, M., Lovasz, L. and Schrijver, A., Eds.) (1988). Geometric algorithms and combinatorial optimization. Springer-Verlag. Volume 2 of Algorithms and Combinatorics. Isidori, A. and C.I. Byrnes (1989). Steady state response, separation principle and the output regulation of nonlinear systems. In: Proc. 28th IEEE Conf. Decision Contr.. IEEE. Tampa, FL. pp . 2247-2251. Mit (1991). ACSL (Advanced Continuous Simulation Language) Reference Manual. Edition 10.0. Natke, H.G. (1988) . Einfiihrung in Theorie und Praxis der Zeitreihen- und Modalanalyse. 2nd ed .. Vieweg. Braunschweig/Wiesbaden. Nemirovsky, A. and D. Yudin (1983) . Problem complexity and method efficiency in optimization . Wiley. Nesterov, Y. and A. Nemirovsky (1993). Interior point polynomial methods in convex programming: Theory and applications. SIAM. Nikolaou, M. (1994) . How nonlinear is a "nonlinear" system? Old and new results under a unifying theory. preprint. Scherer, C. (1994) . Multiobjective H2/ Hoc control. IEEE Trans. Automat. Contr. accepted. Sourlas, D . and V. Manousiouthakis (1992) . Development of linear models for nonlinear plants. AIChE Annual Meeting, Miami, FL. Vidyasagar , M. (1978). Nonlinear Systems Analysis. Prentice Hall. Englewood Cliffs, NJ.
Definition and Computation of a Nonlinearity Measure Frank Allgower* *Inlttitut Jir SYlttemdyno.milc und Rege/unglttechnilc, Univerltitiit Stuttgo.rt, 70550 Stuttgo.rt, FRG, phone: +49-711 - 6856193; Jo.x: +49-711-6856371; emo.il: o.lIgower~rult . uni-lttuttgo.rt.de
Abstract. Nonlinearity measures are a means to quantify the size of nonlinearity in the 1/0behavior of nonlinear systems. This paper is concerned with the definition and the computation of nonlinearity measures. We review some previous results on the computation of nonlinearity measures and also demonstrate a novel approach to efficiently compute nonlinearity measures based on convex optimization. An example illustrates the theory and demonstrates the complexity of the computations involved. Key Words. nonlinearity measures; quantification of nonlinearities; convex optimization
1. Introduction
discrete time systems that can be represented by second order Volterra series. In (Nikolaou, 1994) a different measure is proposed that quantifies the average nonlinearity of a nonlinear I/O-system in discrete time. A completely different nonlinearity measure related to the achievable closed loop performance is described in (Allgower, 1994b). A comprehensive discussion of different measures can be found in (Allgower , 1994a). The practicality of the different definitions depends of course on the application and goal at hand , but also on the computational complexity of its calculation . In this paper the focus is on the derivation of applicable schemes for the computation of a particular nonlinearity measure that appears to be of practical usefulness. We will review our previous results on the calculation of this nonlinearity measure and introduce a new method for its computation based on convex optimization. The paper is structured as follows: In Section 2 the nonlinearity measure is introduced and discussed. Section 3 is concerned with the computation of this nonlinearity measure . In Section 3.1 calculation of a lower bound based on the computation of Fourier coefficients is described . Section 3.2 describes a novel approach based on convex optimization. A short discussion about possible applications of nonlinearity measures in connection with nonlinear control system analysis and design together with an example is given in Section 4.
Not all nonlinear systems have the same "degree" of nonlinearity. There are systems with nonlinearities that have no significant influence on the I/O-behavior. Other nonlinear systems have a behavior that is very different from the behavior of any linear system. By definition , both types of systems are nonlinear. However the first type of systems could be termed "mildly" nonlinear, whereas the second type could be called "strongly" nonlinear. For analysis and synthesis purposes it is often crucial to distinguish between them. In many cases a distinction simply based on intuition is not sufficient. The introduction of nonlinearity measures is an attempt to systematically approach this problem by quantifying the size of nonlinearity in the I/O-behavior of nonlinear systems. The term nonlinearity measure was first used in (Desoer and Wang , 1981) . Without the need for an explicit computation, it was shown there that certain linear feedback laws have a linearizing effect on nonlinear systems. The first systematic application of nonlinearity measures to the analysis of nonlinear control systems and the synthesis of nonlinear controllers was reported in (Allgower and Gilles, 1992) together with algorithms for their computation . Independently nonlinearity measures were also applied in (Nikolaou, 1994) and (Sourlas and Manousiouthakis, 1992) where also computational schemes for their calculation are developed. The measure considered in (Sourlas and Manousiouthakis , 1992) is equivalent to the measure defined in (Desoer and Wang, 1981) and to the one used in this paper. The computational scheme presented in (Sourlas and Manousiouthakis, 1992) is however restricted to
2. Definition of a nonlinearity measure We consider general nonlinear (i.e. not necessarily linear) , multivariable I/O-systems with input u and output y . The input signals u are
257
assumed to belong to the normed space U a of admissible input signals. The output signals y are elements of the normed space of output signals y. The I/O-behavior is described by the nonlinear operator N that maps input signals u E Ua into output signals y E Y :
Other sensible definitions for non linearity measures can be defined as well. For example a normalization of (3) by the "size" of y (= N[u)) may be useful in some cases (Allgower, 1994b): ~u N
y = N[u).
(1)
lim yet)
=0
for
u(t)
==
0
"It.
Definition 1 The nonlinearity measure
sup
In general the computation of nonlinearity measure
(2)
GEy UEU
IIC[u]- N[u]llpy , lIull puG
(3)
3.1. A lower bound on
The "size" of the nonlinearity of an I/O-system is thus defined as the normalized largest difference between the output of nonlinear system Nand a linear system C measured by the norm 11 . Ilpy in y. The difference is taken with respect to the worst-case input in U and with respect to the linear system C* that approximates N best. An example for a specific choice of nonlinearity measure is given by restricting U to consist of L 2 -signals with bounded energy lIullL 2 :S u and the norms 11 . IIpy and 11 . IIpu G both being the L 2-norm . The set U contains the subset of inputs that we consider relevant in the controlled loop. It is easy to see that the nonlinearity measure
Restriction ofU to a suitable set can greatly facilitate the computation of the nonlinearity measure. For reasons of clearer presentation we restrict ourselves in this subsection to SISO systems with continuous inputs and outputs u, y E C[O, 00) satisfying
with C : Ua ........ Y being a linear I/O-operator, 1I·lIpy a norm in y , 11·llpua a norm in U a , and U being the set of inputs considered.
IIC*[u]- N[u]lIpy = 0
(5)
3. Computation of the nonlinearity measure
of a BIBO-stable, causal I/O-system N : U a ........ Y for inputs signals u E U ~ Ua is defined by the nonnegative number ",U . f 'l'N= 10
IIG[u) - N[u)llpy IIN[u]lIpy .
In the remainder of this paper we restrict our attention to the nonlinearity measure defined by (3).
We consider static as well as dynamic operators N, but restrict our attention to causal, BIBOstable systems in this paper. This setup describes a very general class of systems that comprise finite dimensional systems in continuous and discrete time as well as infinite dimensional systems. Without loss of generality, we assume t-+oo
. f
= J~y ;~ft
limT-+oo lIullvT limT-+oo lIyllvT
< <
00
(6)
00
with
lIyllVT =
liT
-
T
y2(t) dt.
(7)
0
We consider the special class of sine-inputs u(t) = A . sin(wt),
(8)
where the amplitude A and frequency ware taken from the sets
A EA = {A E R+ 10:S A :S Amar} (9) wEn = {wER+ IWmin:SW:SWmar}. This specific set of inputs, that we will denote by U$ in the sequel, is thus parametrized by A and w :
(4)
U.
= {u 1 u(t) = A ·sin(wt),A E A,w En} .
(10)
For linear systems the dynamic behavior is completely described by the response to sine-inputs with one arbitrary amplitude and w E R+. Despite the additional variation of the amplitude, the behavior of a nonlinear system N is of course not completely described by the response to input signals in U•. Any nonlinearity measure based on
The behavior of nonlinear system N for the input signals considered is thus the same as that of the linear system C* and therefore N behaves linear in a certain sense. We will refer to this as: "N is linear in U". It must be stressed that the value of the nonlinearity measure depends of course on the choice of input signals U considered .
258
U. instead of U will only be a lower bound on the measure we are really interested in. In particular we suggest the following approximate nonlinearity measure X~' for N:
r
IIG[u] - N[u]lIvT lIullvT
If we consider the set of input signals U. instead of the single input (8), then by the parametrization of U. the relation to be proved follows immediately. 0
(15)
Calculation of the lower bound (11) thus only involves computation of Fourier coefficients for different periodic output signals. For simple nonlinearities, like for example static memory less nonlinearities, an analytical calculation is possible in analogy to the methods for determining describing functions (e.g. Gelb and Vander Velde (1968)). For systems that are more complex, the Fourier coefficients, and thus the lower bound, can be calculated in an efficient and reliable way by numerical techniques (e.g. Natke (1988)). The numerical computation of X~' is only based on the simulated output signal. Thus no restrictions have to be assumed on the type of nonlinear systems that can be considered. Therefore for any system that can be simulated, including systems described by differential-algebraic equations, systems with nonsmooth nonlinearities , etc. , this lower bound on the nonlinearity measure can be calculated. From (11) it follows immediately that if N is a linear I/O-operator, the value of the lower bound is also equal to zero. From X~' = 0 it can however only be concluded that for each input signal u E U there exists a linear system Gu (depending on the respective u), so that the outputs of Nand Gu do not differ for this specific input signal u.
Theorem 1 If the set U. is given by (1 0) and the measure X~' is defined as in eq. (11), then
3.2. Numerical computation of the nonlinearity measure through convex optimization
U, A
XN =
.
f
:~E, J~gT:'rr;.,
(11)
Note that sup and inf are interchanged in (11) as compared to (3). This, together with U. ~ U, results in X~' being a lower bound on 4>#
< A-.U _'I'N,
U,
XN
where
1I ·lIp»
and
1I ·llpu
a
in
(12)
4># are also defined as (13)
The advantage of this particular lower bound is that it can be calculated in a very efficient way. Under some technical assumptions, the response of nonlinear system N to a sine-input (8) with fixed A E A and fixed wEn is the sum of two functions (Isidori and Byrnes, 1989): (i) the periodic steady-state response YIO , and (ii) the transient response Ytr that decays to zero, because of the assumed stability of N : YN
= N[u] =
Ytr
+ YIO ·
(14)
The periodic signal YIO can thus be expanded in a Fourier series with amplitude coefficients Ak 00
YIO
= Ao +
2: A
k'
sin(kwt
+ 4>10) .
10=1
U
XN'
= sup AEA wEn
1 -.
A
2( A) Ao w,
~ AHw , A)
+~
2
In this section we want to demonstrate an alternative method to approximately compute the nonlinearity measure 4>}j. defined by (3) . This approach is based on solving an appropriate nonlinear programming problem numerically. In order to apply numerical methods we have to make two approximations:
' (16)
10=2
where Ak(W , A) are the amplitude coefficients of the ph harmonics of the steady-state response YlO to inputs u(t) = A . sin(wt) .
sketch of proof: As a first step we consider only one input signal ofform (8) with fixed amplitude A and fixed frequency w . From the principle of harmonic balance (e.g. Vidyasagar (1978)) we know that the output YL of the linear system G* that approximates N best in the sense of eq.(l1) is given by the first harmonic of the output of N yL(t)
= G*[u](t) = AI ' sin(wt + 4>1)'
• We restrict the set of inputs considered to a finite set, that we will call Uc in the sequel. • We predefine a fixed structure for the unknown best linear approximation that only contains a finite number of variable real parameters. We will make no further restriction on Uc neither with respect to the type of inputs considered (sinefunctions, steps, random inputs, etc.) nor with respect to the number of input signals. The problem we want to solve is: Find G* E g, such that
(17)
with AI , 4>1 as in (15). Using this G* , it follows after some straightforward algebraic manipulations that inf lim IIG[u] - N[u]llvT GEgT-oo
(18)
= lim IIG*[u] - N[u]lIvT = T-oo
< 259
where 9 is the space of linear operators mapping u E Ua into y .
very satisfying results. Further theoretical investigations on the appropriate choice of basis functions and with respect to convergence properties are however needed. Calculation of the approximate solution B~c of the nonlinearity measure ,p via numerical optimization as described has some desired properties due to the convexity of the problem (Nemirovsky and Yudin, 1983):
Theorem 2 The infinite dimensional optimization problem (19) is a convex programming problem. proof: We have to show that
f :9
~
#-
R+ with
f(G) = max IIG[u]- N[u]lIpy UEUc lI u ll pU4
(20)
is a convex function with respect to G and that 9 is a convex set . As 9 is a Banach space, convexity of 9 follows immediately. In order to show the convexity offunction f we define another function fft with ii E Ue fixed and
fft(G)
= IIG[ii] -
N[ii]llpy·
(21)
This function is convex with respect to G by the definition of convexity: 11
A Gdii] + (1-A)G 2[ii]- N[ii]lIpy (22) = IIAGdii] + (1-A)G2[ii]- AN[ii] - (l-A)N[ii]lIpy ::; AIIG1[ii] - N[ii]lIpy + (1-A)IIG 2 [ii] - N[ii]lIpy-
Because the maximum of a finite number of convex functions (scaled by constants lIullpU4) is again a convex function (Ekeland and Temam, 1976), f(G) = maxftEUc lIulfpU4 . fft(G) is also a convex function. 0
Therefore 8~c can be calculated numerically in an effective and reliable way. If the set of inputs Ue considered for calculation of B~c is not dense in U , then B~c will only be a lower bound on ,p#-. Further investigations will concern the use of an appropriate choice of basis functions for Ue , such that U e is dense in U. The computation of 8~c is more involved than computation of the lower bound X~·. If only sineinputs u E Ue = U. (compare (10» are considered and the norm in (11) and (19) is equally defined for both cases , then X~· ::; B~· . Thus B~' is a tighter bound on ,p~' than X~· . Another advantage of the proposed numerical solution via convex optimization as compared to calculation of X~· is that we are not restricted to consider only sineinputs. Any input signals that appear to be appropriate for the problem at hand can be used for the computation of B~c.
The convex infinite dimensional optimization problem (19) can be approximated through a convex finite dimensional problem by a suitable convex parametrization of the space of linear operators. A possible parametrization for discrete time SISO systems is for example given by assuming a moving average structure (step response model) for linear system G : n
G[u](k) =
LXi ·u(k - i) ,
(23)
;=0
where the scalars Xi are the step response coefficients to be determined . This choice preserves the convexity between the variable vector of parameters Z = [Xl , .. . , xn]T and function fin (20). For continuous time SISO systems for instance a convex parametrization
4. Example
n
G[u] =
LXi ·G;[u],
• There is a unique global optimum. No further local optima exist. • The cost to determine this global optimum increases the most polynomially with the number of variable parameters. • The cost grows the most linearly with respect to the number of input signals u contained in Ue . • In order to determine the global optimum, the response of nonlinear system N to each input u E U e has to be computed only once. Thus a calculation for complex systems is as feasible as a calculation for a simple system. • There exist very powerful algorithms to solve convex optimization problems (e.g. (Nesterov and Nemirovsky, 1993; Grotschel et al. , 1988» .
(24)
Possible applications of nonlinearity measures for the analysis of nonlinear control systems are of course manifold and include control structure selection, validation of nonlinear models and calculation of the best linear approximation. Besides their importance for a general analysis of nonlinear systems, nonlinearity measures are especially interesting in connection with approximate 1/0linearization where they can be used as a tool for a quantitative assessment of the size of nonlinear-
;=1
with G;[u] fixed stable linear systems , is also possible. The extension to multivariable systems is straightforward. It is clear that the quality of the approximate solution depends on the choice of approximation functions. Our own experience and experience collected in a different context (Boyd et al. , 1988; Scherer , 1994) shows however that such finite dimensional parametrizations lead to
260
To demonstrate the flexibility of the approach, we base the computation of eJ:je on a discretized model of (25) that is derived using finite forward differences with relatively fast sampling 6T = 0.2. Considering a set of inputs of 24 different stepfunctions with amplitudes between -0.6 and +0.6
ity that remains after approximate linearization. In this context, nonlinearity measures can also be used to determine the number of approximation terms needed in order to achieve a satisfying linearity. For a detailed discussion see (Allgower, 1995; Allgower, 19946). In this section we want to quickly compare the computational methods presented using a simple example. This will also serve to demonstrate the numerical effort needed and to discuss implications of the choice of input signals U . We consider the following nonlinear SISO system with unstable zero dynamics Y=
Uc1
for the computation of eJ:je leads to a value of
(30) for the uncompensated system. h(t) in (29) denotes the unit step function . The convex optimization was carried out using the MINIMAX function of the MATLAB Optimization Toolbox (Grace, 1992). For this system, with Uc1 as in (29) , a simulation horizon of 100 timesteps, and a step response model with n = 25 coefficients, about ten seconds CPU time were needed on the same VAX-Station as above. It is also possible to compute e in continuous time. For this we consider 24 pulse functions of different amplitude
(25)
with scalar input u, scalar output y, and state :I: = [Xl , X2JT . We are interested in the behavior of
:- ········:··.:.:11' ::.::.:.·:.::.::.:..:::.: .. ::::·:1[·.·:':".::"'::.::','··············· .....·1j I
:
i
i
...................
I
,
I
,
i.................,,.
et
i
.. _-_ ... ,;
!
··· .... ·: .. ·.. ········· .. ·t· .. ·· .. ··· .. ·· .. ·, ··· ···· .. ·...... \
!
with
A E {-0.6, -0.55, .. . , +0.55, +0.6} } (29)
Xl
Xl = -3XI + 4X2 + 2xr + u X2 = -Xl + X2 - 0.5x~,
= {ulu(t) = A . h(t)
u = 0.2
.. ...... !
u
~--+--+--- i
u
= 0.4
= 0.6
Uc 2
= {ulu(t) = A . h(t) -
A . h(t - 50) with
A E {-0.6, -0 .55 , ... , +0.55, +0.6}} Fig. 1. Step responses of uncompensated systems (25) .
and use the parametrization of the best linear approximation according to (24) with only five fixed, stable linear systems
system (25) in a neighborhood around the steady state :1:$ = [O,OV and assume that the input is constrained by
lul:::; 0.6
C;[u]
(26)
.
xt·
eJ:je2 = 0.67.
(33)
Computation time was about four seconds. Using convex optimization it is also possible to calculate 1 with respect to the input signals described by (27). The resulting measure has a value of
et·
= {u I u(t) = A· sin(wt) with A E {0.3, 0.6} and w E {10, l,O.5,O.l,O .05 , 0.01} }
X·
r. ' ,T; E {.5,.99,1.73,2.82,4.45} . (32) s ; +1
Optimization was carried out on the same VAX computer as above using ACSL and an (almost) standard SQP algorithm, leading to a value for the nonlinearity measure of
System (25) displays strongly nonlinear 1/0behavior as demonstrated by the step-responses shown in Figure 1. 1 Calculation of the lower bound of the nonlinearity measure for N using a set of twelve sine-inputs U. I with different amplitudes and frequencIes U. I
(31)
(27)
eJ:j·1 = 0.62
leads to a value of
(34)
for the continuous as well as for the discrete setup . The resulting best linear approximations C* of N derived by all three methods, virtually show identicall/O-behavior. From these figures it can be seen that the methods to compute the nonlinearity measure presented in Section 3.1 and Section 3.2 lead to very similar values of the nonlinearity measure. The 1 value of the lower bound only differs from eJ:jel and eJ:je2 by less than ten percent, sustaining our experience that the lower bound is usually rather tight. It can also be seen from the above that the nonlinearity measure is rather insensitive
(28) for system (25) . The value of the nonlinearity measure, i.e. the gain from input u to the difference between the output of nonlinear system N and the output of its best linear approximation, is thus almost as big as the gain of the nonlinear system itself. This confirms our observation that system N is strongly nonlinear. Computation of 1 for system (25) using the set of inputs (27) in the simulation environment ACSL (Mit, 1991) only needed fractions of a second on a VAX-Station 4090.
xt·
xt·
261
to the choice of input signals as long as the relevant operating regime and dynamic range is captured properly. Even the type of parametrization for the best linear approximation doesn't seem to alter the results in a significant way. System (25) can be approximately 1/0linearized by feedback u = -1.9xr
+ 0.05x~ -
0.14xIX2
+v
(35)
despite its unstable zero dynamics (Allgower, 1995). With this feedback the non linearity measure of the compensated system is
OJj·l =
X~;l = 0.05 .
(36)
Figure 2 shows step responses of the compensated system (25) with feedback (35). The considerably "more linear" behaviour of the compensated system is clearly reflected by the value of the nonlinearity measure. .......... : ........ .
....... :
. ..... .1
..... 1. ......... ..i ...............L · i !
j
= =
v -0 .6 v = -0.4 v -0 .2
v = 0.2 v = 0.4 v = 0.6
Fig. 2. Step responses of system (25) with compensation (35) .
5. Conclusions In this paper we discussed a measure for quantification of the size of nonlinearity in the 1/0behavior of a large class of dynamical systems. This nonlinearity measure captures the difference in the I/O-behavior between the nonlinear system and its best linear approximation. Computation of this nonlinearity measure is usually rather involved. We have given an overview over two methods to approximately calculate the nonlinearity measure: a lower bound based on evaluation of Fourier coefficients, and a novel approach using convex optimization . These approaches are computationally very attractive. Even for rather complex systems of high order a numerical solution is possible. Acknowledgement: Fruitful discussions with C. Scherer are greatfully acknowledged. We also want to thank D. Odenthal, W . Schobbert, and M. Klander for their invaluable help in implementing the computational codes.
6. REFERENCES Allgower, F . (1994a). Entwurfvon naherungsweise E/ A-linearisierenden Reglern fur nichtlineare Systeme. Interner Bericht. ISR, Universitat Stuttgart. Allgower, F. (1994b). NichtlinearitatsmaBe - Ein Werkzeug zur Analyse und Synthese ni(:htlinearer Regelkreise. In : Entwurf Nichtlinearer
262
Regelungen (S. Engell, Ed.). Oldenbourg Verlag. Allgower, F. (1995). Definition and computation of a non linearity measure and application to approximate I/O-linearization. Internal Report 95-1. ISR, Universitat Stuttgart. Allgower, F. and E.D. Gilles (1992). Approximate input/output-linearization of nonlinear systems. AIChE Annual Meeting, Miami, FL. Boyd, S., V. Balakrishnan, C. Barratt , N. Khraishi, X. Li, D. Meyer and S. Norman (1988). A new CAD method and associated architectures for linear controllers. IEEE Trans. Automat. Contr. AC-33(3), 268-283. Desoer, C.A. and Y.-T. Wang (1981). Foundations of feedback theory for nonlinear dynamical systems. IEEE Trans. Circ. Syst. CAS27(2), 104-123. Ekeland, I. and R. Temam (1976). Convex analysis and variational problems. North-Holland . Amsterdam . Gelb, A. and W . E. Vander Velde (1968) . Multiinput Describing Functions and Nonlinear System Design. McGraw-Hill. New York, NY. Grace, A. (1992). Optimization Toolbox for use with MATLAB. The MathWorks, Inc .. South Natick, MA. Grotschel, M., Lovasz, L. and Schrijver, A., Eds.) (1988). Geometric algorithms and combinatorial optimization. Springer-Verlag. Volume 2 of Algorithms and Combinatorics. Isidori, A. and C.I. Byrnes (1989). Steady state response, separation principle and the output regulation of nonlinear systems. In: Proc. 28th IEEE Conf. Decision Contr.. IEEE. Tampa, FL. pp . 2247-2251. Mit (1991). ACSL (Advanced Continuous Simulation Language) Reference Manual. Edition 10.0. Natke, H.G. (1988) . Einfiihrung in Theorie und Praxis der Zeitreihen- und Modalanalyse. 2nd ed .. Vieweg. Braunschweig/Wiesbaden. Nemirovsky, A. and D. Yudin (1983) . Problem complexity and method efficiency in optimization . Wiley. Nesterov, Y. and A. Nemirovsky (1993). Interior point polynomial methods in convex programming: Theory and applications. SIAM. Nikolaou, M. (1994) . How nonlinear is a "nonlinear" system? Old and new results under a unifying theory. preprint. Scherer, C. (1994) . Multiobjective H2/ Hoc control. IEEE Trans. Automat. Contr. accepted. Sourlas, D . and V. Manousiouthakis (1992) . Development of linear models for nonlinear plants. AIChE Annual Meeting, Miami, FL. Vidyasagar , M. (1978). Nonlinear Systems Analysis. Prentice Hall. Englewood Cliffs, NJ.