A method for controller parameter estimation based

0 downloads 0 Views 782KB Size Report
5. Inman DJ (2006) Vibration with control. Wiley, Chichester. 6. Balchen JG, Andresen T, Foss BA (2003) Reguleringsteknikk,. 5th ed. Department of Engineering ...
A method for controller parameter estimation based on perturbations

Magne Bratland, Bjørn Haugen & Terje Rølvåg

Engineering with Computers An International Journal for SimulationBased Engineering ISSN 0177-0667 Engineering with Computers DOI 10.1007/s00366-013-0312-3

1 23

Your article is protected by copyright and all rights are held exclusively by SpringerVerlag London. This e-offprint is for personal use only and shall not be self-archived in electronic repositories. If you wish to selfarchive your work, please use the accepted author’s version for posting to your own website or your institution’s repository. You may further deposit the accepted author’s version on a funder’s repository at a funder’s request, provided it is not made publicly available until 12 months after publication.

1 23

Author's personal copy Engineering with Computers DOI 10.1007/s00366-013-0312-3

ORIGINAL ARTICLE

A method for controller parameter estimation based on perturbations Magne Bratland • Bjørn Haugen • Terje Rølva˚g

Received: 24 August 2011 / Accepted: 3 January 2013  Springer-Verlag London 2013

Abstract Simulation and prediction of eigenfrequencies and mode shapes for active flexible multibody systems is an important task in disciplines such as robotics and aerospace engineering. A challenge is to accurately include both controller effects and flexible body dynamics in a multidisciplinary system model appropriate for modal analysis. A method for performing modal analyses of such systems in a finite element environment was recently developed by the authors. On issue is, however, that for engineers working in a finite element environment, the controller properties are not always explicitly available prior to modal analyses. The authors encountered this problem when working with the design of a particular offshore windmill. The controller for the windmill was delivered in the form of a dynamic link library (dll) from a third party provider, and when performing virtual testing of the windmill design, it was of great importance to use the ‘‘real’’ controller in the form of the provided dll, rather than re-model it in for instance Simulink or EASY5. This paper presents a method for estimating the controller parameters of PID-type controllers when solving the closed-loop eigenvalue problem for active flexible multibody systems in a finite element environment. The method is based on applying incremental changes, perturbations, to relevant system variables while recording reactions from other system variables. In this work, the theory of the method is derived and the method is tested through several numerical examples.

M. Bratland (&)  B. Haugen  T. Rølva˚g Department of Engineering Design and Materials, Norwegian University of Science and Technology, Richard Birkelands veg 2B, 7491 Trondheim, Norway e-mail: [email protected]

Keywords Modal analysis  Finite element method  Control system  Parameter estimation  Perturbation

1 Introduction Modal analysis and dynamic simulation of active flexible multibody systems—from now on referred to as active mechanisms—are a multidisciplinary challenge. The dynamic performance of such products is strongly dependent on an optimal interaction between the controllers and the mechanical components. An important tool in the optimization of such products is modal analysis, which predicts modal parameters, i.e. natural frequencies, mode shapes and damping ratios, for the active system. Due to the complexity of the mechanical components, both in form and in function, it may be practical to handle such systems through a finite element (FE) approach. Effective time domain dynamic simulations of multibody systems in an FE environment have been described by, for instance, Ge´radin and Cardona [1] and Sivertsen [2]. The authors have recently developed a method for performing modal analyses of active mechanisms in an FE environment [3]. In that work, the equations for the control system are expressed in second-order form, rather than in first-order or state-space form, which is typical practice in control system disciplines; see for instance [4–6]. One of the advantages of this approach is an increased compatibility with the mechanical equations, which are typically expressed in second-order form, e.g. [2, 7–10], since equations determined in state-space form are difficult to transform into second-order structural dynamics equations [11]. The generalized eigenvalue problem will be of size n, where n is the number of degrees of freedom (DOFs) and traditional FE eigenvalue problem solvers can be utilized.

123

Author's personal copy Engineering with Computers

One remark about that method is that the controller properties have to be known explicitly prior to the modal analysis. To the best of the authors’ knowledge, no commercial software system for simulation of active mechanisms fully integrates flexible multibody dynamics and control system simulation, since the equations for control systems are typically expressed in first-order form ðx_ ¼ Ax þ Bu; y ¼ Cx þ DuÞ while they are for mechanical systems typically expressed in second-order form ðM€r þ C_r þ Kr ¼ FÞ: Control system software, such as MATLAB and Simulink,1 usually support both controller design and control system simulation where the mechanical system can be modeled with rigid bodies, lumped masses, inertias, springs, dampers or analytical equations. This will cause the flexible body dynamics to be predicted by very simplified models. In flexible multibody dynamics software systems, such as FEDEM,2 feedback type controllers will typically calculate loads applied to the mechanism based on feedback measurements of the system [2]. In addition, some flexible multibody dynamics software systems also have the option of importing or communicating with the controller model as an external process, for instance, through a dynamic link library (dll) or Simulink. For these reasons, the controller is comparable to a ‘‘black box’’ or unknown function, as seen from the mechanical part of the software system. This approach works well in a time domain analysis when the controller drives the mechanism with applied loads based on the given controller algorithms, however, a major problem occurs in modal analyses of the closed-loop system. In free vibration analysis, all loads are set to zero, which decouples the controller and mechanical model. As a result, the mechanism becomes singular in all controlled DOFs. In order to overcome this issue, methods for identifying the controller parameters may be applied. This paper is focused on presenting a method for estimating controller parameters for systems containing either higher-order integral gains, higher-order derivative gains or a combination of proportional, integral and derivative gains, the latter often being referred to as a proportional-integralderivative (PID) controller, the most common type of controllers in use today [13, 14]. However, the method presented in this work is not limited to apply to such controllers only. It may also be applied to any system containing properties corresponding to the ones listed above. An example can be the identification of the 1

MATLAB and Simulink by The MathWorks, Inc. FEDEM (Finite Element in Dynamics of Elastic Mechanisms) simulation software is a multibody dynamics package distributed by Fedem Technology AS. It is based on the finite element method and uses model reduction techniques to effectively perform nonlinear time domain dynamic simulations of active flexible multibody systems [2, 12].

2

123

properties of a mechanical system equal to mass, damping and stiffness based on, for instance, position, velocity or acceleration parameters only. The objective of this work is to derive a method which can be used when performing modal analyses of active mechanisms, using FE based software systems. The software systems used in this paper are MATLAB and Simulink3 and FEDEM.4

2 Interaction between mechanism and controller The equation of motion for a single degree of freedom (SDOF) mechanical system with a single-input singleoutput (SISO) feedback controller can be written as [3]: _ þ k rðtÞ ¼ FApp ðtÞ þ FCtrl ðtÞ m r€ðtÞ þ c rðtÞ

ð1Þ

where m is the mass, c is the damping and k is the stiffness. r is the displacement of the mass m with respect to time; r_ and r€ are the first and second time derivatives of r, i.e. velocity and acceleration of the mass m. FApp is the applied mechanical force and FCtrl is the force from the controller. This is in accordance with equations found in [15]. Figure 1 shows a simple block diagram used for describing a SISO feedback control system. In Fig. 1, y0 is the reference variable, y is the measured variable and e is the difference between y0 and y. u is the controller output and FCtrl is a force from the controller exerted by an actuator. x is the state variable from the physical process (i.e. position r, velocity r_ or acceleration r€), and v is the disturbance on the physical process. Only feedback controllers will be dealt with in this work, hence all control system terminology used here refers implicitly to feedback controllers. For a feedback PID-type controller, the controller output u is given by: Z d uPID ðtÞ ¼ Kp eðtÞ þ Ki ð2Þ eðtÞ dt þ Kd eðtÞ dt where Kp is the proportional gain, Kd is the derivative gain and Ki is the integral gain from the controller. Since e is the difference between y0 and y, the controller output can be split into a feedforward or feedthrough part governed by y0 and a feedback part governed by y, as shown in [16]. The feedforward part can be interpreted as an applied force whose parameters are not affected by the system itself and will not affect the internal dynamics of the system. Therefore, it is not of particular interest in this context. The only part which does affect the internal dynamics of the system is the feedback part. Thus, Eq. (2) can more conveniently be written as: 3 4

MATLAB and Simulink version R2010a. FEDEM version R5.0.

Author's personal copy Engineering with Computers Fig. 1 Block diagram for a SISO feedback control system

v y0

e

Controller

u

Actuator

FCtrl

Physical Process

x

Sensor

y



uPIDFeedback ðtÞ ¼ Kp yðtÞ þ Ki

Z

yðtÞ dt þ Kd

d yðtÞ dt

ð3Þ

One view of the control system is to isolate the control elements from the physical process. The control elements then principally contain three parts: a sensor, an actuator and a controller that contains the various controller elements, as shown in Fig. 2. As shown in Fig. 2, the effects by the control elements on the mechanical system can be given as: oFCtrl oFCtrl ou oy ¼ ox ou oy ox

or dFCtrl ¼ GAct GCtrl GSens dx ð4Þ

where GAct is the actuator gradient, GCtrl is the controller gradient and GSens is the sensor gradient. Similarly, the gradients for a multiple-input multipleoutput (MIMO) system can be written as: dFCtrli ¼

oFCtrli ouj oyk dxl ¼ GActij GCtrljk GSenskl dxl ouj oyk oxl

ð5Þ

or, in matrix form, as: dFCtrl ¼

oFCtrl ou oy dx ¼ GAct GCtrl GSens dx ou oy ox

ð6Þ

Hence, as explained in [3], the equation of motion for the free vibration of a multiple degree of freedom (MDOF) mechanical system with a MIMO feedback controller can thus be written as: M€rðtÞ þ C_rðtÞ þ KrðtÞ þ GAct GCtrl GSens xðtÞ ¼ 0

ð7Þ

where M is the n  n mass matrix, C is the n  n damping matrix, K is the n  n stiffness matrix and r; r_ and €r are the n  1 position, velocity and acceleration vectors, respectively. x is a vector of the system state variables, that is, position, velocity and acceleration. GAct ; GCtrl and GSens are the actuator gradient, controller gradient and sensor gradient matrices, respectively. The actuator gradient GAct describes the relationship between the controller forces FCtrl exerted by the actuator and the output signals u from the controller, and has

dimensions nFCtrl  nu where nFCtrl is the number of controller forces and nu is the number of controller outputs. The controller gradient GCtrl describes the relationship between the input variables y and output variables u both to and from the controller, respectively; that is, the various controller gains. Matrix GCtrl has the dimensions nu  ny where nu is the number of controller outputs and ny is the number of controller inputs. The sensor gradient GSens describes the relationship between the controller input variables y and the system state variables r; r_ and €r represented by the vector x, and has dimensions ny  3nr , where ny is the number of controller inputs and nr is the number of all system DOFs. x is given as: 2 3 r x ¼ 4 r_ 5 ð8Þ €r Vector x has the dimensions 3nr  1 where nr is the number of all system DOFs. Each sensor is limited to measure only one state variable in only one single system DOF or between two system DOFs. The matrix product G of the gradient matrices GAct ; GCtrl and GSens has dimensions nFCtrl  3nr : If G is pre-multiplied with the topology matrix relating each controller force FCtrli with its respective system DOFs, and then split into 3 nr  nr matrices, GPos ; GVel and GAcc ; one for each state variable r; r_ and €r; the matrices GPos ; GVel and GAcc can be added to their respective system matrix yielding the following equation system for the free vibration of a controlled mechanism [3]: ðM þ GAcc Þ€rðtÞ þ ðC þ GVel Þ_rðtÞ þ ðK þ GPos ÞrðtÞ ¼ 0 ð9Þ As shown in Eqs. (7) and (9), there is an uncoupling between the mechanical system and the control elements, hence, the properties of the mechanical system are not of relevance when concerned with identifying the controller parameters.

3 Estimation of controller parameters Control elements y0 (t ) x(t )

Sensor

y (t )

Fig. 2 Control elements

Controller

u (t )

Actuator

FCtrl (t )

One of the main motivations behind this paper is to be able to perform accurate modal analyses of active mechanisms using FE based software systems. To be able to do so, the various controller gains, and hence the controller’s equivalent mechanical properties, must be known. However,

123

Author's personal copy Engineering with Computers

these values are not always explicitly available for the person performing the modal analysis since mechanical engineers and control engineers usually operate in different software systems. Consequently, a method for estimating the values of interest should be derived. A potential method for parameter estimation is to introduce perturbations into the system. This approach is not to be confused with the perturbation method described in [7], which can be used to solve nonlinear differential equations in which the solution is in the form of a power series. Perturbations in this context are incremental changes in a system variable. The basis of this technique can be found in, for instance, the principle of virtual work [7, 10, 17], the displacement method/direct stiffness method [18], system identification/parameter estimation [19, 20] and optimization theory [21]. For all the various fields listed above, the concept remains the same: apply changes in one variable, measure reactions from other variables and then, process the results to derive the desired system parameters. In this paper, perturbations will be used on the decoupled controller to estimate the desired controller parameters. The controller is treated as a ‘‘black box’’ or an unknown function. By applying incremental changes, perturbations, to the input of the controller, small changes in the output from the controller can be registered. These changes will be in accordance with the internal control routine of the controller. The parameters of the controller can thus be estimated based on predetermined changes in the controller input and registered changes from the controller output. One important feature of the proposed technique is a save-and-restore capability of the system variables. After perturbing, all system variables are reset to their pre-perturbation state to not affect any other simulations.

In Fig. 3, the variables time t and controller input y are perturbed by the values Myj and Mtj during perturbation j. y0 and t0 are the initial values for y and t, respectively, at the present time step. From Fig. 3, the following relationships can be derived: Myj ¼ yj  y0 ) yj ¼ y0 þ Myj

ð10Þ

Mtj ¼ tj  t0 ) tj ¼ t0 þ Mtj

ð11Þ

Since the controller output u is a function of y and t, the following equation can be given: Muj ¼ uj  u0 ¼ uj ðyj ; tj Þ  uðy0 ; t0 Þ

The values Myj and Mtj can be chosen arbitrarily, but it can be practical to express Myj as a function of Mtj : The linear equation for yj ðtÞ for perturbation j can then be written as: Myj t ð13Þ yj ðtÞ ¼ y0 þ Mtj 3.2 Estimation of controller parameters for controllers containing proportional gain For a feedback type controller containing only a controller output u proportional to the input variable y, its feedback gain equation can be written as: uP ðtÞ ¼ Kp yðtÞ

A perturbation, as described in the previous section, is illustrated in Fig. 3. y (t )

ou dy or du ¼ Kp dy oy

ou My or Mu ¼ Kp My oy

ð16Þ

Kp can thus be calculated by solving the following equation:

Δy j

Kp ¼ ðMyÞ1 Mu y0 t t0 Fig. 3 Perturbation j of t and y

123

ð15Þ

Since the perturbation technique is meant to be used with computers, it is more suitable to treat Eq. (15) numerically rather than analytically. In discrete differential form, Eq. (15) can be written as: Mu ¼

yj

ð14Þ

where Kp is the proportional gain. Equation (14) can be written on a general differential form as: du ¼

3.1 The perturbation technique

ð12Þ

¼ uj ðy0 þ Myj ; t0 þ Mtj Þ  uðy0 ; t0 Þ

Δt j

tj

ð17Þ

A perturbation algorithm for estimating Kp can be broken into six steps. Since the controller parameters may not be constant with time, the perturbation algorithm should be performed each time an eigenvalue analysis is to be performed. The steps in the perturbation algorithm are:

Author's personal copy Engineering with Computers

1. 2.

Obtain the initial values y0 and u0 for the controller. Establish Myj and Mtj . For simplicity, Myj can be given as Myj ¼ Mtj ;, making it sufficient to establish Mtj . Calculate yj and tj in accordance with Eqs. (10) and (11). Iterate the controller with these new values for the input yj and time tj ; and record the reaction from the controller uj due to the change in the input. Calculate Muj based on u0 and uj in accordance with Eq. (12). Use Eq. (17) to estimate Kp.

3. 4.

5. 6.

uII ðtÞ ¼ Kii

3.3.1 Single integration For a feedback type controller containing only a controller output u proportional to the time integral of the input variable y, its feedback gain equation can be written as: Z ð18Þ uI ðtÞ ¼ Ki yðtÞ dt where Ki is the integral gain. Equation (18) can be written in discrete differential form as: Z Z ou M y dt or Mu ¼ Ki M y dt Mu ¼ R ð19Þ o y dt and Ki can be calculated by solving the equation:  Z 1 Muj Ki ¼ M yj dt

ð20Þ

In order to estimate Ki using the perturbation technique R described in the previous sections, M yj dt needs to be R discretized. In Fig. 3, M yj dt is the area under the linear R curve. If yj ðtÞ is given as in Eq. (13), M yj dt can be made as a function of Myj and Mtj by: Z M

  1 Myj 2 Mtj yj dt ¼ y0 t þ t 2 Mtj 0 0   1 ¼ y0 þ Myj Mtj 2

yj dt ¼

ZMtj

Inserting Eq. (21) into Eq. (20) yields:   1 1 Ki ¼ y0 þ Myj Mtj Muj 2

ð21Þ

ð22Þ

yðtÞ dt dt

ð23Þ

where Kii is the double integral gain. Equation (23) can be written in discrete differential form as: ZZ ZZ ou M y dt dt or Mu ¼ Kii M Mu ¼ RR y dt dt ð24Þ o y dt dt M

Based on antidifferentiation, the double integral RR yj dt dt can be derived in discrete form as: ZZ

M 3.3 Estimation of controller parameters for controllers containing integral gain

ZZ

  1 2 1 Myj 3 Mtj yj dt dt ¼ yj dt dt ¼ y0 t þ t 2 6 Mtj 0 0 0   1 1 y0 þ Myj Mtj2 ¼ ð25Þ 2 6 ZMtj Z t

Kii can thus be estimated by solving:   1 1 1 2 y0 þ Myj Mtj Muj Kii ¼ 2 6

ð26Þ

3.3.3 Triple integration Like the single and double integral presented in the previous sections, the feedback gain equation for a feedback type controller containing only a controller output u proportional to the triple time integral of the input variable y can be written as: ZZZ uIII ðtÞ ¼ Kiii yðtÞ dt dt dt ð27Þ where Kiii is the triple integral gain. Following the same procedure as for the double integral in the previous section, the RRR triple integral M yj dt dt dt can be derived in discrete form as: ZZZ ZMtj Z t Z t M yj dt dt dt ¼ yj dt dt dt 0 0 0  1 3 1 Myj 4 Mtj t ¼ y0 t þ 24 Mtj 0 6 1 1 3 y0 þ Myj Mtj ¼ 6 24 and Kiii can be estimated by solving:   1 1 1 y0 þ Myj Mtj3 Muj Kiii ¼ 6 24

ð28Þ

ð29Þ

3.4 Estimation of controller parameters for controllers containing derivative gain

3.3.2 Double integration 3.4.1 Single derivation Like the single integration presented in Sect. 3.3.1, the feedback gain equation for a feedback type controller containing only a controller output u proportional to the double time integral of the input variable y can be written as:

For a feedback type controller containing only a controller output u proportional to the time derivative of the input variable y, its feedback gain equation can be written as:

123

Author's personal copy Engineering with Computers

uD ðtÞ ¼ Kd

d _ yðtÞ ¼ Kd yðtÞ dt

ð30Þ

where Kd is the derivative gain. Equation (30) can be written in discrete differential form as: Mu ¼

ou My_ or Mu ¼ Kd My_ oy_

ð31Þ

For a feedback type controller containing only a controller output u proportional to the double time derivative of the input variable y, its feedback gain equation can be written as: uDD ðtÞ ¼ Kdd

Kd can be calculated by solving the equation:  1 Kd ¼ My_j Muj

ð32Þ

In order to estimate Kd using the perturbation technique described in the previous sections, My_ j has to be discretized. The derivative in Fig. 3 can be given as: Myj y_ j ¼ Mtj

ð33Þ

d2 yðtÞ ¼ Kdd y€ðtÞ dt2

ð34Þ

However, y_ 0 does not exist in Fig. 3. In order to have both My_ j and My_ 0 , two perturbation steps have to be performed. An example of a two-step perturbation is illustrated in Fig. 4. As for the one-step perturbation illustrated in Fig. 3, the values My0 ; Myj ; Mt0 and Mtj can be chosen arbitrarily, but it can be practical to express Myj as a function of Mtj ; while Mt0 can be given as Mt0 ¼ Mtj and My0 ¼ 0. Equation (34) can then be simplified to: Myj Mtj

ð35Þ

Mu ¼

ou M€ y or Mu ¼ Kdd M€ y o€ y

ð36Þ

ð38Þ

Similarly to the discretization of My_ j in Sect. 3.4.1, M€ yj can be written as: My_j My_ 0  Mtj Mt0

ð39Þ

In order to derive M€ y; three perturbation steps have to be performed. An example of a three-step perturbation is given in Fig. 5. In Fig. 5, Myj0 is the first perturbation step, Myj1 is the second perturbation step and Myj2 is the third perturbation step in perturbation j. In the figure, the following parameters are given: Mtj0 ¼ Mtj1 ¼ Mtj2 ¼ Mtj ; Myj0 ¼ 0; Myj1 ¼ Mtj and Myj2 ¼ Myj1 : Using this three-step perturbation series with the parameters as shown in Fig. 5, yj can be given as: the variables Myj ; My_j and M€ Myj ¼ Myj2 My_j ¼ My_ j2 ¼

ð40Þ Myj2  Myj1 Mtj My My

Inserting Eq. (35) into Eq. (32) yields:  1 Myj Kd ¼ Muj Mtj

ð37Þ

where Kdd is the double derivative gain. Equation (37) can be written in discrete differential form as:

M€ yj ¼ y€j  y€0 ¼

Using Eq. (10) as a basis, My_j can be given as: Myj My0 My_j ¼ y_ j  y_ 0 ¼  Mtj Mt0

My_j ¼ y_ j  0 ¼

3.4.2 Double derivation

ð41Þ

j2 j1  My_ j2  My_ j1 Mtj M€ yj ¼ ¼ Mtj Mtj Myj2  2 Myj1 þ Myj0 ¼ Mtj2

Myj1 Myj0 Mtj

ð42Þ

y (t )

y (t ) yj

y j1

Δy j

Δy j1

Δy j 2

Δy j 0

y0 Δy0 y−1

y j2

t

t−1 Δt0 Fig. 4 Two-step perturbation

123

t0

Δt j

tj

Δt j 0

Δt j1

Fig. 5 Three-step perturbation

Δt j 2

t

Author's personal copy Engineering with Computers

Since Myj0 ¼ 0; Eq. (42) can be reduced to: M€ yj ¼

Myj2  2 Myj1 Mtj2

Kdd can thus be calculated by solving: !1 Myj2  2 Myj1 Kdd ¼ Muj Mtj2

ð43Þ

To avoid singularities when performing the matrix inversion in Eq. (49), the determinant of the invertible matrix should be nonzero. This requirement is met for My1 6¼ My2 6¼ My3 and Mt1 6¼ Mt2 6¼ Mt3 . Typically, Mtj and Myj can be given as:

ð44Þ

Mt1 ¼ d  Mtsim ; Myj ¼ j  My1

3.5 Estimation of controller parameters for controllers containing combinations of proportional, integral and derivative gains For a PID controller, the feedback controller output is given by Eq. (3). In a discrete differential form, this can be written as: Z ou ou ou R Mu ¼ My þ M y dt þ My_ or Mu oy oy_ o yZdt ¼ Kp My þ Ki M y dt þ Kd My_ ð45Þ As demonstrated in Eq. (45), a PID controller is a compound controller consisting of both a proportional gain, an integral gain and a derivative gain. In order to estimate all three gains Kp, Ki and Kd using the perturbation technique, three perturbations need to be performed. And since the controller contains a derivative gain, a two-step perturbation algorithm needs to be used, as explained in Sect. 3.4.1. This gives the following set of equations: Z Mu1 ¼ Kp My1 þ Ki M y1 dt þ Kd My_ 1 Z Mu2 ¼ Kp My2 þ Ki M y2 dt þ Kd My_ 2 ð46Þ Z Mu3 ¼ Kp My3 þ Ki M y3 dt þ Kd My_ 3 which can be written in matrix form as: R 2 3 2 32 3 Mu1 My1 M R y1 dt My_1 Kp 4 Mu2 5 ¼ 4 My2 M y2 dt My_2 54 Ki 5 R Mu3 My3 M y3 dt My_3 Kd

Mtj ¼ j  Mt1 ;

ð50Þ

where Mtsim is the simulation time increment. d is a small positive scalar called the relative perturbation step size [21]. A possible default value of d; as used by the authors in this work, is 0.1. A perturbation algorithm for estimating Kp, Ki and Kd can be broken down into eight steps. Since the controller parameters may not be constant with time, as mentioned in Sect. 3.2, the perturbation algorithm should be performed each time an eigenvalue analysis is to be performed. The steps in the perturbation algorithm are: 1. 2. 3. 4. 5. 6.

7. 8.

Do one initial perturbation on the controller with My ¼ 0 and Mt 6¼ 0: This is to ensure y_ 0 ¼ 0: Obtain the initial values y0 and u0 for the controller. Establish Myj and Mtj . For a PID controller, j ¼ 1. . .3; but if Eq. (50) is used, it is sufficient to establish Mt1 : R Calculate M yj dt and My_ j . These are given by Eqs. (21) and (35). Calculate yj and tj in accordance with Eqs. (10) and (11). Iterate the controller with these new values for the input yj and time tj ; and record the reaction from the controller uj due to the change in the input. Calculate Muj based on u0 and uj in accordance with Eq. (12). Use the matrix system in Eq. (49) to estimate Kp, Ki and Kd.

3.6 Partitioning of perturbation steps ð47Þ

To derive the controller properties Kp, Ki and Kd, one can solve the following matrix system by the use of matrix inversion: R 2 3 2 31 2 3 Kp Mu1 My1 M R y1 dt My_ 1 4 Ki 5 ¼ 4 My2 M y2 dt My_ 2 5 4 Mu2 5 ð48Þ R Kd My3 M y3 dt My_ 3 Mu3 Inserting Eqs. 2 3 2 My1 Kp 4 Ki 5 ¼ 6 My 4 2 Kd My3

My1 ¼ Mt1 ;

(21) and (35) into Eq. (48) yields: 31 2   3 1 y þ 1 My Mt My Mu1 Mt1  0 21 1  1 My 7 y þ My Mt Mt22 5 4 Mu2 5 ð49Þ  0 21 2  2 My Mu3 y0 þ 2 My3 Mt3 Mt33

When testing the perturbation technique presented in Sects. 3.1–3.5, the authors experienced a problem with some controller simulation algorithms. In some cases, the controllers were suffering from what seemed as an erroneous time step dependent delay from when a change in the input resulted in a change in the outputs. In order to extract the gradients from such a system, it may be necessary to perform a perturbation using multiple time steps, that is, introducing incremental steps in between each perturbation. Such a multistep perturbation is illustrated in Fig. 6. As illustrated in Fig. 6, if the perturbation is divided into n equal sections, each perturbation step Mtj and Myj would be composed of n incremental sub-steps Mtji and Myji , where i ¼ 1. . .n. Note that the subscripts in this section are

123

Author's personal copy Engineering with Computers

y (t )

t j3 , y j3

yj

Fig. 7 System of triple integration modeled in Simulink using the Continuous Integrator block

t j2 , y j2

Δy j

y0

Δy j2

t j1 , y j1 t j0 , y j0

Δt j2

Table 1 Comparison of the perturbation technique against results from the simulation in Simulink for a system with triple integration Perturbation

Simulink

Mt1

0.1

0.1

My1 R M y1 dt RR M y1 dt dt RRR M y1 dt dt dt

0.1

0.1

0.005

0.005

t

t0

Δt j

tj

Fig. 6 Partition of perturbation

not to be confused with the subscripting presented in Sect. 3.4.2. For each tji and yji ; the system would be iterated. For the perturbation technique itself, only tjn and yjn would be used. The incremental sub-steps would then be of size: Mtji ¼

1 Mtj n

ð51Þ

4 Testing of the perturbation technique In order to verify the theory and methods derived in Chapter 3, some numerical tests were performed on three different controllers. The first controller was one containing a higher-order integral gain. The second controller contained a higher-order derivative gain, while the last controller was a compound controller containing combinations of proportional, integral and derivative gains (PID-type controller). For the first two cases, the objective was to test the derived theory by comparing it against a commercial software system, represented here by Simulink. Since the perturbation technique is meant to be implemented in FEDEM, which is a FE bases software system, the final case was focused on the accuracy of parameter estimation for variants of PID controllers. 4.1 Testing of the perturbation technique for parameter estimation of controllers containing integral gain A comparison was made between the perturbation technique and Simulink for a system containing single, double and triple integration. For the perturbation technique, Eqs. (21), (25) and (28) were used to calculate R RR RRR M y dt; M y dt dt and M y dt dt dt; respectively. In Simulink, this system was created using three integration blocks in series, as shown in Fig. 7. For the triple integration system, the parameters Mtj ; Myj and y0 were given as Mt1 ¼ 0:1 s, My1 ¼ Mt1 and y0 ¼ 0. In Simulink, the ode4 (Runge–Kutta) solver was used. The

123

1:6667  104

1:6667  104

6

4:1667  106

4:1667  10

simulation start time was set to 0.0 and the simulation stop time to 0.1, with a fixed-step size of 0.1. A comparison between the perturbation technique and the simulation results from Simulink is shown in Table 1. As can be seen in Table 1, the perturbation technique R and the Simulink simulation are identical for M y1 dt; RR RRR M y1 dt dt and M y1 dt dt dt. 4.2 Testing of the perturbation technique for parameter estimation of controllers containing derivative gain A comparison was made between the perturbation technique and Simulink for a system containing single and double derivation. For the perturbation technique, Eqs. (41) and (42) were used for My_ and M€ y; respectively. In Simulink, this system was created using two derivation blocks in series, as shown in Fig. 8. For the double derivation system, the parameters of interest used in the simulations were Mtj0 ¼ Mtj1 ¼ Mtj2 ¼ Mtj ; Myj0 ¼ 0; Myj1 ¼ Mtj and Myj2 ¼ Myj1 , with the value for Mtj given as Mt1 ¼ 0:1 s: In Simulink, the ode4 (Runge–Kutta) solver was used. The simulation start time was set to 0.0 and the simulation stop time to 0.2, with a fixed-step size of 0.1. A comparison between the perturbation technique and the simulation results from Simulink is shown in Table 2. As can be seen in Table 2, the perturbation technique and the Simulink simulation are identical for both My_ 1 and M€ y1 : 4.3 Testing of the perturbation technique for parameter estimation of PID controllers To verify the theory and method derived in Sect. 3.5, some basic tests were performed using the perturbation technique. The objective of the tests was to verify whether the

Author's personal copy Engineering with Computers

2

Fig. 8 System of double derivation modeled in Simulink using the Continuous Derivative block

Table 2 Comparison of the perturbation technique against results from the simulation in Simulink for a system with double derivation Perturbation Mt1

0.1

Simulink 0.1

My1

-0.1

My_ 1

-2

-0.1 -2.0

M€ y1

-30

-30.0

perturbation technique could be used to estimate the controller parameters for any PID-type controller during any time step of a nonlinear dynamic time domain simulation. Since the perturbation technique is intended to be implemented in FEDEM, it is vital to test and verify the method in this software system. In addition, one initial test of the perturbation technique for PID controllers was performed in Simulink. One possible setup for a PID controller in Simulink is shown in Fig. 9. For the PID system in Simulink, the following controller parameters were used: Kp = 1, Ki = 1 and Kd = 1. These are the values which are to be treated as the unknown parameters of the PID controller; although to verify the perturbation technique, the values in this example are known a priori. For the perturbation of the PID controller based on Fig. 4, the following parameters were used: Mtj ; Myj ; Mt0 ; My0 ; t0 and y0 were given as Mt1 ¼ 0:1 s, Mtj ¼ j  Mt1 ; Myj ¼ Mtj ; Mt0 ¼ Mt1 and My0 ¼ t0 ¼ y0 ¼ 0. The ode4 (Runge–Kutta) solver was used, and the simulation start time was set to 0.0. Three perturbations were performed (j = 1…3), and for each perturbation j, the fixed-step size was set to Mtj and the simulation stop time to 2  Mtj . Based on Eq. (49), the results from the simulation in Simulink are shown in Eq. (52).

Kp

3

2

0:1000

0:0050

6 7 6 4 Ki 5 ¼ 4 0:2000 0:0200 0:3000 0:0450 Kd 2 3 1:0000 6 7 ¼ 4 1:0000 5 1:0000

1:0000

31 2

1:1050

3

7 6 7 1:0000 5 4 1:2200 5 1:0000 1:3450 ð52Þ

As revealed in Eq. (52), Kp = 1, Ki = 1 and Kd = 1. These estimated values for Kp, Ki and Kd are identical to the actual values for the controller gains of the PID controller, indicating a validity of the perturbation technique for such controllers. Since the initial simulation results from Simulink indicate that the perturbation technique is valid for PID controllers, three additional tests were performed using the perturbation technique in FEDEM. As previously mentioned, the objective of these tests was to establish whether the perturbation technique could be used to estimate the controller parameters for any PID-type controller during any time step of a nonlinear dynamic time domain simulation. The setup for the tests is shown in Fig. 10. The setup in Fig. 10 consists of a SDOF system with mass m, damping c and stiffness k. There is only one DOF: position r of the mass. r is the input for the controller, which is of type PID. As stated in Sect. 2, the properties of the mechanical system are not of relevance when concerned with identifying the controller parameters, hence, neither the parameters value nor the number of DOFs of the mass-spring-damper system are of relevance in this context. The numerical values for the controller gains were arbitrarily chosen, but were deliberately given different values to more easily distinguish between them. The values for Mtj and Myj were given by Eq. (50), with d ¼ 0:1: The simulation time increment Mtsim for all tests in this section was set to 0.01 s, i.e. Mt1 ¼ 0:001: Three different tests were performed on the active massspring-damper system shown in Fig. 10. The first test was to insure that the perturbation technique worked for controllers of any possible combination of P, I and D. The next test was to insure that the perturbation technique worked at any time step of the dynamic time domain simulation and not only at the start-up of the simulation. The last test was to insure that the perturbation technique would also work for discontinuous systems. 4.3.1 Various controller combinations of P, I and D

Fig. 9 PID controller modeled in Simulink

The perturbation technique should be able to yield correct estimations of the controller parameters for any PID-type controller for the possible combinations of P, I and D, including the trivial system without any controller present. To verify this, the active system in Fig. 10 was set in static

123

Author's personal copy Engineering with Computers Fig. 10 SDOF Mass-springdamper system with position feedback PID controller

r

Kp

k

v

y0

m

c

FCtrl

e

K i ∫ dt

u

Actuator

FCtrl

mr + cr + kr

r

Sensor

y

− Kd

d dt

Table 3 Estimated controller parameters for different combinations of P, I and D controllers Kp

Ki

Kd

PID

100.000000000005

19.999999999069

6.000000000000

PI

100.000000000000

19.999999999884

0.000000000000

PD

100.000000000000

0.000000000000

6.000000000000

ID

0.000000000000

20.000000000000

6.000000000000

P

100.000000000000

0.000000000000

0.000000000000

I

0.000000000000

20.000000000000

0.000000000000

D None

0.000000000000 0.000000000000

0.000000000000 0.000000000000

6.000000000000 0.000000000000

equilibrium, and the perturbation technique was performed on a total of eight different types of PID controllers: PID, PI, PD, ID, P, I, D and zero-gain controller. The basic values for the different gains were set to Kp ¼ 100; Ki ¼ 20 and Kd ¼ 6; in addition to zero when not included. The results from the different perturbations are shown in Table 3. The numerical values of the estimations are given here with 12 decimals. With regard to the intended usage of the perturbation technique, correct values up to the second decimal should be more than sufficient, but to depict the accuracy of the method and show when the estimated values deviate from the correct ones, 12 decimals were used. The results presented in Table 3 demonstrate that the perturbation technique yields approximately correct controller parameter estimations for any PID-type controller. All but the PI and PID controller yield correct values up to the 12th decimal. For the PI and PID controller, the integral gain Ki is only correct up to the 9th decimal, and the proportional gain Kp for the PID controller is incorrect on the 12th decimal. 4.3.2 Perturbations on PID controller during time simulation with sinusoidal input signal The perturbation technique should be able to yield correct estimations of the controller parameters for any value of the initial values y0 and u0 for the controller, i.e. yield correct estimations of the controller parameters at any time step of the dynamic time domain simulation. To verify this, the position r of the mass m in the active system in Fig. 10 was given a prescribed sinusoidal motion of 1 Hz, and the perturbation technique was performed on a PID controller at

123

Fig. 11 Input signal to the controller. The crosses mark each time the perturbation sequence is performed

various time steps. The simulation ran for one second with a time increment of 0.01 s. The perturbation technique was performed both at start-up and at time intervals of 0.1 s, giving a total of 11 different perturbations. The gains to the PID controller were set to Kp ¼ 100; Ki ¼ 20 and Kd ¼ 6. The input signal to the controller is shown in Fig. 11. The results from the simulation are shown in Table 4. As can be seen in Table 4, the perturbation technique yields correct results up to the 8th decimal at any time step for any of the controller parameters for this simulation. Both the Kp and Ki estimations are accurate up to about the 8th decimal, while the Kd estimations are accurate up to the 13th decimal. For the proportional gain, Kp, the perturbation technique yields the greatest errors at time 0.2 and 0.7 s, both being approximately 2:0  108 : For the integral gain, Ki, the perturbation technique yields the greatest errors at time 0.1, 0.2 and 0.7 s, the difference being approximately 2:5  108 for 0.1 s and approximately 2:0  108 for 0.2 and 0.7 s, respectively. For the derivative gain, Kd, the perturbation technique yields the greatest errors at time 0.1 and 0.2 s, being about 5  1014 and 4  1014 ; respectively.

Author's personal copy Engineering with Computers Table 4 Estimated controller parameters at different time steps for sinusoidal input signal Ki

Time

Kp

0.0

100.000000000005 19.9999999990687

0.1 0.2

99.999999986962

Kd 6.00000000000000

20.0000000223517

6.00000000000005

100.000000018626 19.9999999795109

5.99999999999996

0.3

99.999999991618

20.0000000083819

0.4

99.999999997206

20.0000000037253

6.00000000000003

0.5

100.000000000000 20.0000000018626

6.00000000000000

0.6 0.7

99.999999994412 19.9999999916181 100.000000019558 20.0000000204891

5.99999999999999 6.00000000000003

0.8

99.999999994412

6.00000000000003

19.9999999944121

6.00000000000002

0.9

100.000000000000 20.0000000009313

5.99999999999999

1.0

100.000000000005 19.9999999990687

6.00000000000000

4.3.3 Perturbations on PID controller during time simulation with discontinuous sinusoidal input signal To further test the capabilities of the perturbation technique to yield correct estimations of the controller parameters for any value of the initial values y0 and u0 for the controller, the system described in Sect. 4.3.2 was used with a discontinuous input signal. The sinusoidal signal was given a switch-to-zero-value at ±0.7, meaning it would go to zero whenever the absolute value of the input signal became larger than 0.7. As shown in Fig. 12, this should occur at time 0.2, 0.3, 0.7 and 0.8 s. All other parameters were the same as they were in Sect. 4.3.2. The results from the simulation are shown in Table 5. As can be seen in Table 5, the perturbation technique yields correct results up to the 6th decimal at any time step for any of the controller parameters for this simulation. Both the Kp and Ki estimations are accurate up to about the 6th decimal, while the Kd estimations are accurate up to the 12th decimal. For the Kp estimation, the greatest error is encountered at 0.3 s, being approximately 8  107 . For the Ki estimation, the greatest error is encountered at 0.5 s, being approximately 2  106 . For the Kd estimation, the greatest errors are encountered at 0.2, 0.4 and 0.5 s, all being approximately 2:5  1012 .

5 Discussion The main motivation behind this work is to make engineers working in an FE environment able to perform accurate modal analyses of active mechanisms. Today, FEDEM has the capability of performing accurate time domain simulations using the controllers to drive the FE model with applied loads based on the given controller algorithms. These controller algorithms can be created

Fig. 12 Input signals to the controller. The crosses mark each time the perturbation sequence is performed. The input signal is a discontinuous sinusoidal signal with switch-to-zero-value at ±0.7

Table 5 Estimated controller parameters at different time steps for discontinuous sinusoidal input signal Time

Kp

Ki

Kd

0.0

100.00000000000400

19.99999999813740

0.1

100.00000011222400

19.99999981001020

6.00000000000060

0.2

-0.00000000044409

0.00000000000000

0.00000000000222

0.3

-0.00000082446987

0.00000086508578

0.00000000000040

0.4

99.99999936856320

20.00000107102100

6.00000000000230

0.5

99.99999999535250

20.00000215135510

6.00000000000253

0.6

99.99999972432850

19.99999952781950

5.99999999999841

0.7

-0.00000004406638

-0.00000004640732

-0.00000000000005

6.00000000000000

0.8

-0.00000004406638

-0.00000004640732

-0.00000000000005

0.9

100.00000010058300

20.00000017136340

5.99999999999991

1.0

100.00000000017100

19.99999991711230

5.99999999999981

either in FEDEM’s Control Editor or in an external software system, such as for instance Simulink. Therefore, the controller algorithms are not required to be known for the engineer working in the FE environment. As stated in the introduction, a major obstacle for modal analysis of the closed-loop system is that in free vibration analysis all loads are set to zero, thereby resulting in a decoupling of the controller and mechanical model. By identifying the controller parameters, the controller’s mechanically equivalent properties can be added to the FE model for the modal analyses, thus, including both controller and

123

Author's personal copy Engineering with Computers

mechanical properties and hence, improving the accuracy of the modal analysis. The method derived in this work is intended to be implemented in FEDEM, though as presented here, it is not dependent on any particular software system. The results presented in Sects. 4.1 and 4.2 show that the perturbation technique yields correct estimations for systems containing single, double and triple integration, as well as single and double derivation compared to results derived using Simulink. This indicates validity of the perturbation technique for such systems. Controllers containing either triple integration or double derivation are not a very common type of controllers; however, in this work they do serve the purpose of testing the robustness of the perturbation technique, which only strengthens the validity of the derived method for its intended use. The validity of the perturbation technique for PID controllers was briefly tested in Simulink. The results from that initial test demonstrated that the perturbation technique is able to correctly estimate such controllers in Simulink. More thorough tests of the technique were conducted in FEDEM, and as can be seen from the results presented in Sect. 4.3, the perturbation technique yields estimations for controller parameters with a highly satisfactory accuracy for any PID-type controller during any time step of a nonlinear dynamic time domain simulation in FEDEM. The greatest estimation error in any of the tests still yielded correct resulting up to the sixth decimal. This should be more than sufficient since a requirement for satisfactory accuracy should be correct results up to the second decimal with regard to the intended usage of the perturbation technique. Hence, for PID-type controllers, the derived method should be able to provide accurate estimations of the controller parameters. Still, one note about the perturbation technique should be made. There has to be a correlation between the state variables of the perturbed system and those used in the perturbation technique. For instance, it can be tempting to believe that the perturbation technique is able to accurately predict the effective mass, stiffness and damping values for an active system containing a position feedback PID controller by perturbing the active system and deriving the system parameters only with respect to mass, stiffness and damping, and not including the integral gain from the controller. By ignoring some of the state variables of the perturbed system, critical and vital system information can either be lost or estimated to incorrect values, rendering the technique virtually useless for its intended use. However, when used properly, the technique has the potential of being of great assistance in identifying the unknown parameters of a controller, such as when performing modal analysis of active mechanisms in an FE environment.

123

6 Conclusion In this paper, a method for controller parameter estimation by the use of perturbations has been presented. The theory for perturbation of systems containing single, double or triple integral functions, single or double derivative functions or a combination of proportional, integral and derivative functions has been derived and tested using commercial software systems. The results from the tests reveal that the derived theory works well for all the mentioned controller variants. If used properly, the presented technique, with its capabilities of accurate controller parameter estimation, has the potential of being a powerful tool for engineers who are conducting modal analysis of active flexible multibody systems in a finite element environment. Acknowledgments The authors would like to acknowledge the guidance and assistance of Professor Ole Ivar Sivertsen and Professor Kristian Tønder at the Norwegian University of Science and Technology (NTNU). The authors would also like to acknowledge the financial support from the Research Council of Norway and the other partners in the Lean Product Development (LPD) Project.

References 1. Ge´radin M, Cardona A (2001) Flexible multibody dynamics: a finite element approach. Wiley, Chichester 2. Sivertsen OI (2001) Virtual testing of mechanical systems— theories and techniques. Advances in Engineering 4. Swets and Zeitlinger B.V., Lisse, The Netherlands 3. Bratland M, Haugen B, Rølva˚g T (2011) Modal analysis of active flexible multibody systems. Comput Struct 89(9–10):750–761 4. Preumont A (2002) Vibration control of active structures: an introduction, 2nd edn. Kluwer Academic Publishers, Dordrecht 5. Inman DJ (2006) Vibration with control. Wiley, Chichester 6. Balchen JG, Andresen T, Foss BA (2003) Reguleringsteknikk, 5th ed. Department of Engineering Cybernetics, Norwegian University of Science and Technology, Trondheim, Norway (in Norwegian) 7. Thomson WT, Dahleh MD (1998) Theory of Vibration with Applications, 5th edn. Prentice Hall, Inc., Upper Saddle River 8. Palm WJ (2007) Mechanical Vibration. John Wiley & Sons, Inc., Hoboken 9. Cook RD, Malkus DS, Plesha ME, Witt RJ (2002) Concepts and applications of finite element analysis, 4th edn. Wiley, New York 10. Bathe K-J (1996) Finite Element Procedures. Prentice Hall, Englewood Cliffs 11. Alvin KF, Park KC (1994) Second-order structural identification procedure via state-space-based system identification. AIAA J 32:397–406 12. Sivertsen OI, Waloen AO (1982) Non-linear finite element formulations for dynamic analysis of mechanisms with elastic components. In: Washington, DC, USA. American Society of Mechanical Engineers, ASME, New York, NY, USA, p 7 13. Astrom KJ, Hagglund T (2001) The future of PID control. Control Eng Pract 9:1163–1175 14. Astrom KJ, Hagglund T (2004) Revisiting the Ziegler-Nichols step response method for PID control. J Process Control 14:635–650

Author's personal copy Engineering with Computers 15. Alkhatib R, Golnaraghi MF (2003) Active structural vibration control: a review. Shock Vib Dig 35:367–383 16. Bratland M, Rølva˚g T (2008) Modal analysis of lumped flexible active systems (Part 1). In: Paper presented at the SIMS 2008: the 48th Scandinavian Conference on Simulation and Modeling, Oslo, Norway, 07 October 2008–08 October 2008 17. Park KC, Felippa CA, Ohayon R (2009) The d’AlembertLagrange principal equations and applications to floating flexible systems. Int J Numer Meth Eng 77:1072–1099

18. Felippa CA (2001) A historical outline of matrix structural analysis: a play in three acts. Comput Struct 79:1313–1324 19. Astrom KJ, Eykhoff P (1971) System identification-a survey. Automatica 7:123–162 20. Perreault EJ, Kirsch RF, Acosta AM (1999) Multiple-input, multiple-output system identification for characterization of limb stiffness dynamics. Biol Cybern 80:327–337 21. Trier SD (2001) Design optimization of flexible multibody systems. Doctoral Thesis, Norwegian University of Science and Technology

123