Test Cycle Optimization using Regression Analysis

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Test Cycle Optimization using Regression Analysis.

Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping

av Dejen Meless LiTH-ISY-EX--10/4242--SE

Linköping 2010

Department of Electrical Engineering Linköpings tekniska högskola

Linköping University Linköpings universitet

(2)

(3)

Test Cycle Optimization using Regression Analysis.

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Dejen Meless LiTH-ISY-EX--10/4242--SE

Handledare: Shiva Sander-Tavalley

ABB Sweden, Corporate Research

Patrik Axelsson

ISY, LiTH

Examinator: Erik Wernholt

ISY, LiTH

(4)

(5)

Avdelning, Institution

Division, Department ISY

Department of Electrical Engineering Linköping University

SE-581 83 Linköping, Sweden

Datum Date 2010-01-29 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version

ISBN

—

ISRN

LiTH-ISY-EX--10/4242--SE

Serietitel och serienummer

Title of series, numbering

ISSN

—

Titel

Title

Optimering av testcykel med hjälp av regressionsanalys. Test Cycle Optimization using Regression Analysis.

Författare

Author

Dejen Meless

Sammanfattning

Abstract

Industrial robots make up an important part in today’s industry and are assigned to a range of different tasks. Needless to say, businesses need to rely on their machine park to function as planned, avoiding stops in production due to machine failures. This is where fault detection methods play a very important part. In this thesis a specific fault detection method based on signal analysis will be considered. When testing a robot for fault(s), a specific test cycle (trajectory) is executed in order to be able to compare test data from different test occasions. Furthermore, different test cycles yield different measurements to analyse, which may affect the performance of the analysis. The question posed is: Can we find an optimal test cycle so that the fault is best revealed in the test data? The goal of this thesis is to, using regression analysis, investigate how the presently executed test cycle in a specific diagnosis method relates to the faults that are monitored (in this case a so called friction fault) and decide if a different one should be recommended. The data also includes representations of two disturbances.

The results from the regression show that the variation in the test quantities utilised in the diagnosis method are not explained by neither the friction fault or the test cycle. It showed that the disturbances had too large effect on the test quantities. This made it impossible to recommend a different (optimal) test cycle based on the analysis.

Nyckelord

(6)

(7)

Abstract

Industrial robots make up an important part in today’s industry and are assigned to a range of different tasks. Needless to say, businesses need to rely on their machine park to function as planned, avoiding stops in production due to machine failures. This is where fault detection methods play a very important part. In this thesis a specific fault detection method based on signal analysis will be considered. When testing a robot for fault(s), a specific test cycle (trajectory) is executed in order to be able to compare test data from different test occasions. Furthermore, different test cycles yield different measurements to analyse, which may affect the performance of the analysis. The question posed is: Can we find an optimal test cycle so that the fault is best revealed in the test data? The goal of this thesis is to, using regression analysis, investigate how the presently executed test cycle in a specific diagnosis method relates to the faults that are monitored (in this case a so called friction fault) and decide if a different one should be recommended. The data also includes representations of two disturbances.

The results from the regression show that the variation in the test quantities utilised in the diagnosis method are not explained by neither the friction fault or the test cycle. It showed that the disturbances had too large effect on the test quantities. This made it impossible to recommend a different (optimal) test cycle based on the analysis.

(8)

(9)

Acknowledgements

I would like to thank my project colleagues Hans Andersson, Kari Saarinen and Shiva Sander-Tavalley at ABB for all their help and guidance through my work and for giving me the chance to take part in a rewarding experience. I would also like to thank my mentor Patrik Axelsson and my examiner Erik Wernholt for their valuable feedback. Lastly but definitely not the least I would like to thank all of the thesis workers at ABB Corporate Research in Västerås that I had the chance to become friends with, without whom, crunching Matlab code and statistics literature would have been so much harder.

(10)

(11)

Introduction

Robots – man made humanoid machines once created to perform the most tedious and laborious tasks that no man was willing to preoccupy himself with, but one day reached a point where they were able to take control over their own destiny, and did so in a violent rebellion. That’s how the Czech playwright Karel ˇCapek in 1920 told the story about the robota in his play Rossum’s Universal Robots, where the word found its first use, literally meaning work, labour or serf labour1_.

1.1 Definition of a robot

What defines a robot has however not always been a clear subject, many defini-tions have been proposed through the years. The human-like robots that ˇCapek tells about in his play would today be classified as androids. In this report the robots in focus are industrial six-axis robots. The International Organization for Standardization gave in 1994 this definition of a robot

An automatically controlled, reprogrammable, multipurpose, manipu-lator programmable in three or more axes, which may be either fixed in place or mobile for use in industrial automation applications.

In Merriam-Webster’s online dictionary Robotics is defined as technology deal-ing with the design, construction, and operation of robots in automation.

1.2 Fault Detection in robots

The increasing demands on product quality and cost efficiency in the today’s in-dustry leads to ever growing complexity and automation degree. This calls for more safety and reliability in the technical processes. Thus early fault detection and system diagnosis becomes of great importance and is increasingly becoming

1_{http://en.wikipedia.org/ wiki/Robot#Etymology}

(14)

2 Introduction

Figure 1.1. A six-axis industrial robot.

a common ingredient in modern automatic control [2]. This work aims at in-vestigating and optimizing a test cycle to a diagnosis method used for industrial six-axis robots and was performed in collaboration with ABB Corporate Research in Västerås.

1.3 The diagnosis method

The diagnosis method considered is based on fault detection by regular comparison between so called test data and nominal data. For every test, a specific test cycle is executed and measurements of two signals are performed. The signals are then subtracted with their respective nominal (normal) data, forming a residual for each signal. If any fault is present in the system it should thus appear as the difference (residual) between the collected data and the nominal data. The nominal data are measurements from the first tests performed, usually after a robot’s last reset or when it first entered production.

1.3.1 Test cycle

The stated test cycle in the method is to move a robot axis ±5◦ at full motor speed for approximately three seconds. There are no requirements on the position to move the joint about, but is typically around the point 0◦. A test cycle needs to be executed in an identical way every time and needs to be performed for every monitored axis. In this report, however, we will only be analysing axis 3, see Chapter 2 for more on Robotics.

(15)

1.4 Outline 3

1.3.2 Test data

While executing a test cycle the applied motor torque τ (t) and motor angular velocity ˙ϕ(t) are sampled at 248 Hz. After a test has been performed the measured

signals s = [τ ˙ϕ] are respectively synchronised and subtracted with their nominal

values snom = [τnom ϕ˙nom], resulting in a torque and velocity residual sres =

s − snom = [τres ϕ˙res]. In this investigation all the measurements are generated

by a simulator provided by ABB.

1.3.3 Test quantities

When there are residuals ready, two test quantities are computed on the data. A model-fit measure φf it=  1 − q PN k=1(sres(k) − ˆsres(k)) 2 q PN k=1(sres(k) − ¯sres) 2  · 100 (1.1)

measuring correlation in the residual and a root mean square-measure (RMS-measure). φRM S = s PN k=1s2res(k) N ,s_P_N k=1s2nom(k) N · 100 (1.2)

measuring the energy ratio in the residual. Where ˆsres(k) and ¯sresare the

estima-tion and the mean of sres(k) respectively. The variable ˆsres(k) is the estimation

produced by an autoregressive model (AR(n)) of the residual, of order n = 2. The test quantities are computed for each vector in sres(k) and snom(k) (first using

the torque vector and later the velocity vector). Naturally, the test quantities are expected to become larger than nominal values when a fault occurs.

1.3.4 Diagnosis statement

After a period of time of regular testing, the computed values of the test quantities from each test occasion will form a time series of values. There will thus be four time series, a sequence of the model-fit measure computed on the torque signal and a sequence of the model-fit measure computed on the velocity signal and similarly two sequences of the RMS-measure will be available. On each sequence, tests to detect trends and abrupt steps are performed. Based on these tests the method will yield a statement of the robot’s health.

1.4 Outline

The outline of the report will be such that the first chapter introduces the reader to the problem at hand. The following three chapters will lay the theoretical foundation to the analysis that will be performed in the Chapter 7. Chapter 5 aims to familiarise the reader with the analysed signals and a set of disturbance

(16)

4 Introduction models. Chapter 6 presents the method of analysis, and lastly a concluding part and suggestions for future work will final the thesis in Chapter 8.

1.5 Limitations

The method is aimed to be applied on any axis on a robot but we will only be analysing axis 3 in this report, see Chapter 2 for more on Robotics. Furthermore, the analysis will focus on investigating the method’s test quantities and their variation due to fault(s), and leave the part of diagnosis statement aside.

(17)

Chapter 2

Robotics

In this chapter a presentation of the basics in robot kinematics and rigid-body modelling will be given, based on [10]. In Section 2.3 an example of rigid-body modelling using the Euler-Lagrange equations will be given to show the basics in how dynamic modelling of robots can be done.

2.1 Introduction

Figure 2.1. A six axis industrial robot.

The mechanical structure of a robot manipulator is characterised by a set of rigid links interconnected by joints. Joints can be classified as being either pris-matic, performing translational movements, or revolute, performing rotary move-ments. Every joint provides the structure one degree of freedom (DOF). For a robot to be able to execute arbitrary movements in three dimensions six joints

(18)

6 Robotics are required, three for positioning and three for orientation. If a robot has more degrees of freedom it is from a kinematic point of view said to be redundant [9]. Figure 2.1 shows a robot manipulator with six revolute joints.

2.2 Homogeneous transformations

One of the requirements of a robot, as stated in Section 1.1, is that it should be able to move autonomously. This implies the need of a programmed controller that plans the trajectory for the robot’s end-effector (gripper, tool). To do this one must establish the relationship between the trajectory that is to be followed and how to translate this motion to joint rotations. This relationship is referred to as kinematics.

Let’s assume we are looking at a point p described in the two coordinate sys-tems ox0y0z0 and ox1y1z1 as shown in Figure 2.2. Point p has thus different

representations p1and p0 in respective frame.

p0=x0i0+ y0j0+ z0k0 (2.1) p1=x1i1+ y1j1+ z1k1 (2.2) p0 z0 y0 x0 d1 0 z1 y1 x1 p1 p

Figure 2.2. Rotational and translational displacement between two frames.

Vector p1 can be projected onto ox0y0z0, or the base frame, giving the

rela-tionship between the frames as

x0= p0· i0= p1· i0= x1i1· i0+ y1j1· i0+ z1k1· i0 (2.3a) y0= p0· j0= p1· j0= x1i1· j0+ y1j1· j0+ z1k1· j0 (2.3b) z0= p0· k0= p1· k0= x1i1· k0+ y1j1· k0+ z1k1· k0. (2.3c)

(19)

2.3 Rigid body modelling 7 p0= Rp1=   i1· i0 j1· i0 k1· i0 i1· j0 j1· j0 k1· j0 i1· k0 j1· k0 k1· k0     x1 y1 z1  . (2.4) The matrix R, or R1

0, is called the rotation matrix from system ox1y1z1to system ox0y0z0 and represents the coordinate transformation relating the point p in both

frames as a result of the rotational displacement between the frames.

Let P represent the transformation between two points including the transla-tional displacement as well.

P = p 1 (2.5) For any vector p0 in the base frame we have

p0= R10p1+ d10 (2.6)

which, with the new notation, results in the transformation

P0= H01P1= R1 0 d10 0 1 p1 1 , 0 = [0 0 0]. (2.7) The matrix H1

0, is called a homogeneous transformation from frame 1 to frame

0 and describes the simple transformation between two frames.

2.3 Rigid body modelling

In this section we will through an example show how a dynamic model of a robot can be derived, however, not used as a model in this work. The model is aimed to provide a description of the relationship between the applied joint torque and the structure’s motion. Models play important roles for simulating motions and designing control algorithms facilitating a range of analyses before constructing a robot [9].

A way of deriving the equations of motion is to utilise the Euler-Lagrange equations, also referred to as analytical mechanics. In contrast to the Newtonian vector mechanics, analytical mechanics bases the motion relations on calculating the (scalar) kinetic and potential energy. The Euler-Lagrange equations relates a set of so called generalised coordinates qi, which represent a body’s motional

degrees of freedom to a set of generalised forces τi. The equation to solve is

d dt ∂L ∂ ˙qi −∂L ∂qi = τi (2.8) L = K − V, (2.9)

where L is called the Lagrangian and K and V are the kinetic and potential energy respectively.

(20)

8 Robotics Figure 2.3(b) shows a schematic equivalent of axis 2 shown in Figure 2.1. The arm rotates with angle ϕl, the applied motor torque is τ and we have a damping

torque τf = γmϕ˙m+ γlϕ˙l. Furthermore we have the link’s mass M , a distance l

from the point of rotation, and the rotational inertias of the motor and link, Jm

and Jl.

(a) Axis 1 and 2 shown in Figure 2.1.

Jm Jl τ M l τf ϕl ϕm

(b) Schematic equivalent of axis 2.

Figure 2.3. Schematic of axis 2.

We start calculating the kinetic energy

K = Jmϕ˙ 2 m 2 + Jlϕ˙2l 2 (2.10) = 1 2 Jmϕ˙2m+ Jlϕ˙2m n2 (2.11) where ϕl= − ϕm n . (2.12)

The variable n here represents the gear reduction factor. The potential energy becomes

V = M gl cos ϕl= M gl cos

ϕm

n (2.13)

and we get the Lagrangian

L = 1 2 Jmϕ˙2m+ Jlϕ˙2m n2 −M gl cosϕm n . (2.14)

Computing (2.8) treating ϕmas the generalised coordinate q now results in

d dt Jmϕ˙m+ Jlϕ˙m n2 + 1 nM gl sin ϕm n = τ − τf. (2.15)

(21)

2.3 Rigid body modelling 9 Which is rewritten as Jm+ Jl n2 ¨ ϕm+ γm+ γl n ˙ ϕm+ M gl n sin ϕm n = τ. (2.16)

Further simplifying the differential equation, the equation of motion finally results in A ¨ϕm+ B ˙ϕm+ C sin ϕm n = τ (2.17) where A = Jm+ Jl n2 (2.18) B =γm+ γl n (2.19) C =M gl n . (2.20)

(22)

(23)

Chapter 3

Fault detection and isolation

This chapter aims to give the reader a quick review of the field of fault detection and isolation (FDI) and is mostly based on [5, 7, 8]. It is common to include identification when discussing FDI but is nevertheless left out in this discussion.

3.1 Introduction

In every technical application there’s a need to know or estimate the state of a process or a system under observation in relation to the requirements put on the system. These requirements can for example be that the product needs to hold a certain quality or that the process should not be in a state too close to failure. To do this one needs methods to diagnose the system. Figure 3.1 shows a schematic description of a typical diagnosis situation. A system is in operation either autonomously or under control of an operator. The system is constantly subjected to different unknown inputs such as disturbances and faults. The desire is naturally to have the system under constant normal operation, and the duty of a fault-detection system is to monitor the process and generate an alarm if it approaches a state of failure. The earlier a fault is detected the more time is available for taking proper measures.

FDI generally comprise of three functions:

1. Fault detection. To indicate the occurrence of fault(s). 2. Fault isolation. To localise the occurrence of fault(s). 3. Fault identification. To determine the magnitude of the fault(s).

(24)

12 Fault detection and isolation System Control Disturbances Faults Observations Diagnosis system Diagnosis statement

Figure 3.1. A typical diagnosis situation.

Terminology

Disturbance An unknown and uncontrolled input acting on the system. Failure A fault resulting in a permanent interruption of the system’s

performance.

Fault An unacceptable deviation from the nominal state of a system property.

Residual A signal which is ideally zero in the fault-free case or when non-monitored faults are present.

3.2 Limit checking

Traditionally a widely used method of keeping a system under supervision has been to watch for breached thresholds and tolerances in the process. The concept is simple; when a threshold has been breached an alarm is generated and appropriate measures can then be taken. This method is well suited in situations where the system operates approximately around a steady state. The big advantage of limit checking methods is their simplicity [8].

3.3 Hardware redundancy

Another method of detecting and at the same time reducing the risk of system failure is the use of hardware redundancy, in other words to duplicate or triplicate hardware instances in the system, putting for example two sensors instead of a single one to measure the same process variable. In this way a fault is detected by simply comparing one sensor with its redundant instance, and there is no need to stop the process despite the occurrence of a fault. This method is highly reliable and offers direct localisation of the fault in the system. There is however no possibility to determine which sensor of the two that is defective. That can however be overcome by a triplication of hardware. Triple redundancy is common practice in security critical components such as in inner-loop controls of an aircraft [8].

(25)

3.4 Model-based fault detection 13

3.4 Model-based fault detection

The approach of model-based fault detection has been developed since the 1970’s and has developed remarkably since then. Its efficiency in detecting faults has been demonstrated by many successful applications in industrial processes and automatic control systems [2]. In contrast to hardware redundancy where com-parison is made with a physical parallel system, the idea of model-based fault detection, or analytical redundancy, is to use mathematical process models and dependencies of measurable signals to detect faults [8]. Based on measured input signals and output signals, the detection methods generate residuals, parameter estimates or state estimates which are called features. Changes in present and normal features result in analytical symptoms [7]. Model-based fault detection has many advantages by being able to detect smaller faults, facilitating isolation of faults and being a more cost effective alternative to hardware redundancy. A disadvantage may be the need of well designed process models which can pose a difficulty depending on the complexity of the process [5]. Figure 3.2 shows a general schematic of model-based fault detection.

Actuators Process Sensors

Process model Feature generation Change detection u(t) y(t) Features Analytical symptoms Faults Model-based fault detection Nominal behaviour

Figure 3.2. Schematic of general model-based fault detection.

3.4.1 Fault modelling

Modelling fault in a process means in general to try to give an analytical description to the faults that are to be monitored. A common way of modelling faults is to model them as deviations in constant model parameters from their nominal values. The aim of fault modelling is of course to detect a fault before it causes a failure,

(26)

14 Fault detection and isolation so better fault models, where for example smaller faults can be detected earlier, naturally implies better detection performance. Typical faults modelled in this way are “gain-errors” and “biases” in sensors [5]. Faults can in their time behaviour be abrupt, incipient or intermittent faults [7], illustrated in Figure 3.3.

f

t

(a) Abrupt fault. f t (b) Incipient fault. f t (c) Intermittent fault.

Figure 3.3. Time behaviour of faults.

3.4.2 Residuals and residual generation

A fundamental part in model-based diagnosing is the process of formulating resid-uals. A residual can be regarded as a time-varying signal used as a fault detector to reveal the deviation of the system’s current state from its nominal state. As mentioned in Terminology (Section 3.1) the residual is typically designed to be zero when no fault exists and non-zero when a fault does exist. A main difficulty in achieving this is the disturbance decoupling problem, i.e. to obtain a residual free from undesired influences that are not considered as faults [5, 8]. Generally this is very difficult to achieve ideally, so residual generation is typically followed by residual evaluation which can comprise of filtering and testing against limits and tolerances [2].

Residual generators are more or less filters that produce a residual from mea-surable process inputs and outputs aiming to fulfill the above stated requirements. In [7] a few examples of residual design approaches are presented, described in the following sections.

3.4.3 Fault detection using parity equations

A straightforward example of model-based diagnosis is to use a fixed model GM

of the monitored process GP and compute the output error

e(s) = [GP(s) − GM(s)] u(s) (3.1)

to detect any inconsistencies, also referred to as a parity equation. To utilise this method however a priori knowledge of the process is required and for single-input-single-output (SISO) processes, only one residual can be generated and one can not distinguish between different faults [7]. Figure 3.4 shows the scheme of fault detection using parity equations.

(27)

3.4 Model-based fault detection 15 ˆ y(t) u(t) y(t) e(t) − GP GM

Figure 3.4. Fault detection using parity equations (simulation).

3.4.4 Fault detection using signal models

Measured signals can many times contain useful information such as harmonic or stochastic oscillations. Sensors can for example be used to monitor machine vibration to detect imbalances, bearing faults, knocking or chattering. If changes in these oscillations can be related to faults in the process, actuators or sensors, signal analysis can be a further source of information. The Fast Fourier Transform, autocorrelation functions or spectral densities

Xν = N −1 X k=0 xke−i2πν k N R_xx(γ) = g X γ=0 xkxk−γ Φ(ω) = |Xν|2 2π (3.2)

can then be useful tools to analyse data. The spectral density and auto correla-tion funccorrela-tions are in these cases very useful in separating stochastic signals from periodic signals. The deficiencies of these methods, however, are that they are not well suited if the frequency content of the signal is unknown on beforehand. Then it is preferred to use parametric models, which are sensitive to even small frequency changes, in order to estimate the main frequencies and their amplitudes [7].

3.4.5 Observer based fault detection

The basic idea behind the development of the observer-based fault diagnosis tech-nique is to replace the process model by an observer which will deliver reliable estimates of the process outputs as well as to provide the designer with the de-sign freedom to achieve the desired decoupling using the well established observer theory [2]. The method is suitable for multivariable processes as opposed to the parity equations approach [7].

Assume that we have a state space model

˙

x(t) = Ax(t) + Bu(t) (3.3a)

(28)

16 Fault detection and isolation where A, B and C are known and the model is a correct representation of the system. A state observer

˙ˆx(t) = Aˆx + Bu(t) + He(t) (3.4a)

e(t) = y(t) − Cˆx(t) (3.4b)

is used to reconstruct the unmeasurable state variables based on the measured inputs and outputs. The errors are computed as

˙˜

x(t) = ˙x − ˙ˆx (3.5)

˙˜

x(t) = [A − HC]˜x(t). (3.6) The state error ˜x vanishes asymptotically if the observer is stable, which can be achieved by proper design of the observer feedback H [7].

3.5 Change detection

Detecting change is the central task in FDI. But as described in previous sections a process is generally under influence of disturbances and detecting a change in the fault indicators means thus to find a significant change. The measured or estimated quantities such as parameters, state variables or signals are usually stochastic variables Si(t) with mean value and variance

¯

Si= E[Si(t)]; σ¯i2= E[ Si(t) − ¯Si

2

] (3.7)

as normal values in a non-faulty case. Changes, or analytic symptoms, are then obtained as the deviations

∆Si= E[Si(t) − ¯Si]; ∆σi = E[σi− ¯σi] (3.8)

from the normal values. To determine then if a change is deviating significantly enough, methods of hypothesis testing are applied [7, 8].

3.6 Fault isolation

The objective of fault isolation is to determine what type of fault or where a fault has occurred. There are several methods of achieving fault isolation, a few examples from [8] being;

Directional residuals The principle of directional residuals is to design a resid-ual vector that changes direction depending on the fault that is present. One can then by classifying certain faults to specific directions isolate the type of fault that is acting on the system. According to [8] this approach has not been widely used, probably because of the difficulties of designing a residual vector with desired properties.

(29)

3.6 Fault isolation 17 Structured residuals A more used approach is the method of structured resid-uals. Structured residuals are sets of residuals where every residual reacts to a specific subset of faults. By observing which residuals are responding one can in that way isolate the fault.

Structured hypothesis tests In [8] structured hypothesis tests are described as being a generalisation of structured residuals but where the isolation method is more formally defined to be available to work with any fault model.

(30)

(31)

Chapter 4

Modelling of time series

In this chapter a brief presentation of two methods of modelling time series, linear regression and autoregression, will be given. A derivation of the R2_{statistic, F}

-statistic, confidence intervals and the model-fit will be given. The largest part of the chapter is based on [3, 4]. The section on autoregressive models however is based on [6].

4.1 Introduction

Time series are sequences of data points that take on different stochastic be-haviours where the temporal ordering of the data is of crucial importance. Time series are typically measured at equidistant time points and the goal of time se-ries analysis is to find patterns and meaningful statistical properties in the data in question. Methods for time series analysis can be divided in two classes; frequency-domain methods and time-frequency-domain methods. Frequency-frequency-domain methods comprise for example of spectral analysis of time series while time-domain methods comprise mainly of autocorrelation and cross-correlation between data points. A few exam-ples of time series are; share prices on the stock market, temperature fluctuations during a day, or, in our case a residual in a diagnosis process.

4.2 Linear regression

A common tool to model time series is the method of linear regression. Linear regression has come to be used in a wide range of applications and the method basically models the relationship between a measured response variable y and a linear combination of measured regression variables x1, x2, . . . , xn, where y is an

observation of the stochastic variable

Y = β0+ β1x1+ . . . + βnxn+

and is normally distributed with zero mean and variance σ. Furthermore, all k

are independent for every observation k = 1, 2, . . . , N . To estimate the parameters

(32)

20 Modelling of time series

θ = [β0β1. . . βn]Twe use the method of least-squares, which is based on estimating

the parameters so that the squared error between observed y and model output is minimised. In algebraic terms, it should fulfil

min θ N X k=1 (yk− E[Yk])2 (4.1a) E[Y ] = β0+ β1x1+ . . . + βnxn (4.1b)

where E[Y ] stands for the expectation of Y.

Example 4.1: Curve fitting

Let’s assume we have been given the data set and the plot shown in Figure 4.1, where y is a 10×1 vector and x is a 10×2 matrix. The question posed is; Which line

y = ˆβ0+ ˆβ1x1

passes through the data so that there is equal error to data points above the line as well as below the line, and how do we calculate it?

y [x0 x1] 11.61 1 1.00 11.46 1 2.00 10.02 1 3.00 15.68 1 4.00 14.58 1 5.00 y [x0 x1] 17.96 1 6.00 18.28 1 7.00 15.63 1 8.00 15.42 1 9.00 18.85 1 10.00 0 2 4 6 8 10 12 10 12 14 16 18 20 x y

Figure 4.1. Table of paired data and its corresponding plot. x0symbolises the constant that is multiplied with ˆβ0.

We recall from (4.1) Q (β0, β1. . . , βn) = N X k=1 (yk− β0− β1xk1− . . . − βnxkn) 2 . (4.2)

Given n ≤ N we can choose to see this as a system of n equations and consider solving ∂Q

∂ ˆβi

= 0. We thus have a system of equations

∂Q ∂βi = N X k=1 2(yk− β0− β1xk1− . . . − βnxkn)(−xki) = 0 (4.3)

(33)

4.2 Linear regression 21

xT(y − xβ) = 0 (4.4)

or

xTy = xTxβ. (4.5)

These equations are called the normal equations to the solution. If now det(xTx) 6= 0, we get ˆ β = xTx−1xTy (4.6) ˆ β = 10.54 0.80 . (4.7)

4.2.1 The R

2

-statistic

After computing the regression coefficients in the previous example it is useful and also common practice to evaluate the estimated model, performing for example an analysis of variance (ANOVA) and/or computing the R2_{-statistic. The R}2_statistic

indicates how much of the variation in the observed data that can be explained by the model and how much of the variation that is attributed to error. A value R2

= 0 is equivalent to using only β0 to model the response variable, while R2 = 1

corresponds to a perfect model that yields exact predictions. We start deriving

R2 _{by first stating}

yk− ¯y = (yk− ˆyk) + (ˆyk− ¯y). (4.8)

Equation (4.8) simply states that the vertical distance between an observed data point yk and the mean ¯y is the sum of the distances between the point and the

model output ˆyk and between the model output and the mean. It can be shown

that N X k=1 (yk− ¯y) 2 = N X k=1 (yk− ˆyk) 2 + N X k=1 (ˆyk− ¯y) 2 (4.9) or otherwise formulated SST OT = SSE+ SSR. (4.10)

Which says that the total squared sum (SS) of variation is equal to the squared sum of variation related to model error plus the squared sum of variation related to the regression model, and R2 is

R2= SSR

SST OT

(34)

4.2.2 t-statistics

Confidence intervals can be computed for both the response variable and the es-timated regression coefficients for the purpose of estimating their interval of vari-ation and furthermore show a variable’s significance in the model, i.e. show if its interval includes zero (variation around zero indicates insignificance). A confidence interval for β is in an ideal situation, i.e. when the entire population of samples is known, defined as

I = ¯β ± σ.

For non-ideal situations, as in nearly every practical situation, we need to estimate the population mean and variance and make use of Gossets t-distribution. The

t-distribution is the probability function of the ratio A

pB/c where

• A is normally distributed with expectation value 0 and variance 1.

• B is a χ2_{-distribution with c degrees of freedom.}

• A and B are independent, i.e. the probability of observing a certain value in A does not depend on the value of B.

We will derive an expression for the confidence interval of β by first defining a help variable. It can be shown that

ˆ βi∼ N βi, σ p hii (4.12) where xTx−1 =      h00 h01 · · · h0n h10 h11 · · · h1n .. . ... . .. ... hn0 hn1 · · · hnn      .

We can thus state

ˆ

βi− βi

σ√hii

∼ N (0, 1). (4.13)

It can also be shown that

(N − n − 1)s2

σ2 ∼ χ

(35)

4.2 Linear regression 23

where

s2= SSE

N − n − 1.

Furthermore, using the definition of the t-distribution, (4.13) and (4.14), we can state ˆ βi− βi σ√hii , s (N − n − 1)s2_/σ2 (N − n − 1) ∼ t(N − n − 1) (4.15) which can be rewritten as

ˆ

βi− βi

s√hii

∼ t(N − n − 1). (4.16)

Relation (4.16) can now be used to compute the confidence interval I by first looking up the needed t-value T and creating

ˆ βi− T · s p hii≤ βi≤ ˆβi+ T · s p hii. (4.17)

4.2.3 F -statistics

The F -statistic is mainly used to investigate if it is worthwhile to add more pa-rameters to a stated linear model and is based on assessing the magnitude of the extra sum of squares derived from the added parameters. If we assume that the added extra model parameters make no significant difference in the model and that SS_R(1) and SS_R(2) represent the sums of squares from regression for model 1 (original) and model 2 (with added parameters), it can be shown that

w =

SS_R(2)− SS_R(1)/p

SS(2)_R /(N − n − p − 1) (4.18)

is an observation of the stochastic variable

W ∼ F (p, N − n − p − 1),

where p is the number of added parameters. N and n remain as the number of observations and number of original parameters respectively. The larger the sum

SS_R(2) becomes the more likely it is that the added parameters make a significant difference in the model. A common criterion is that w is so large that it belongs to the 10% - 5% largest values in the F -distribution, i.e. so large that there would be 10% to 5% percent chance (or more correctly, risk) of observing a larger w.

(36)

4.3 Autoregressive models

As mentioned previously time series are temporally ordered data sequences that can retain certain characteristics in the time-domain as well as in the frequency-domain. The key to time series modelling is to recognise these patterns, and another way of doing so is to use autoregressive models. Autoregressive models are models that are based on the assumption that present values can be modelled as a combination of past values. The difference equation

y(t) + β1y(t − 1) + . . . + βny(t − n) = α0e(t) + α1e(t − 1) + . . . + αme(t − m) (4.19)

is called an autoregressive moving average-model of order (n,m), also written ARMA(n,m). It is quite common to have the case m = 0, resulting in a so called AR(n)-model [6]

y(t) = −β1y(t − 1) − . . . − βny(t − n) + e(t). (4.20)

Similarly the case of n = 0 results in a so called moving average model MA(m)

y(t) = α0e(t) + α1e(t − 1) + . . . + αme(t − m). (4.21)

To estimate θ = [β1β2 . . . βn] T

in (4.20) we go about in a similar fashion as in the regression example by first stating

E[Y ] = −β1y(t − 1) − . . . − βny(t − n) = ϕ(t)Tθ (4.22)

where

ϕ(t)T = [−y(t − 1) − y(t − 2) . . . − y(t − n)]

θ = [β1β2 . . . βn]T.

Again utilising the least-squares method we state the minimising criterion

min θ 1 N N X t=1 y(t) − ϕ(t)Tθ2 . (4.23)

Skipping the intermediate steps we state that the following can be shown

ˆ θ = " N X t=1 ϕ(t)ϕ(t)T #−1 N X t=1 ϕ(t)y(t). (4.24)

The expression in (4.24) shows that the best approximation ˆθ is based entirely on

(37)

4.3 Autoregressive models 25

4.3.1 Model-fit

A way of assessing autoregressive models is to calculate the model’s goodness of fit in the same way as in Section 1.3.3, and is computed as

V2= 1 − PN t=1(y(t) − ˆy(t)) 2 PN t=1(y(t) − ¯y) 2 , y(t) = ϕˆ T_θ._ˆ _(4.25)

There are a couple of issues to address, mentioned in [6], when estimating autoregressive models.

• Deterministic components such as trends should be checked for before mod-elling a signal to not get misleading results.

• If the signal is oversampled this should also be addressed. A rule of thumb is that the signal should be approximately sampled ten times the bandwidth of the signal.

• A signal may contain high frequency components that are not desired to be modelled. One should consider filtering the signal with an appropriate low-pass filter. It should be noted however that the filter will affect the signal dynamics.

(38)

(39)

Chapter 5

Measurements and

disturbance models

The purpose of this chapter is to familiarise the reader with the robot measure-ments that will be analysed and also describe the disturbances that are represented in the simulated measurements.

5.1 Measurements

Figure 5.1 shows the torque τ (t) and velocity ˙ϕ(t) from a simulation of a test. It

can be seen that the motion is periodic. Figure 5.2 shows the frequency content in these signals, revealing that the signals’ period lies around 2.5 Hz.

0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 Time [s] Torque [Nm]

(a) Torque signal in time domain.

0 0.5 1 1.5 −200 −150 −100 −50 0 50 100 150 200 Time [s]

Angular velocity [rad/s]

(b) Velocity signal in time domain.

Figure 5.1. Measurements in time domain.

(40)

28 Measurements and disturbance models 0 5 10 15 20 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Frequency [Hz] Amplitude

(a) Torque signal in frequency domain.

0 5 10 15 20 0 200 400 600 800 1000 Frequency [Hz] Amplitude

(b) Velocity signal in frequency domain.

Figure 5.2. Measurements in frequency domain.

5.2 Disturbances

As discussed in Chapter 3 a fundamental part of model-based diagnosis is fault modelling. In this report the investigated diagnosis method is evaluated on its ability of detecting a friction fault in a robot joint under the influence of a set of disturbance factors. A description of their representation follows below.

5.2.1 Friction fault

The monitored fault is represented as a changed friction torque and will be referred to as friction fault. The extra torque does not have to be attributed to any specific cause or causes. The investigation is concerned with all faults that imply an added torque of the kind that is explained by the friction model

τf( ˙ϕm) = Fc+ (Fs− Fc) e−( ˙ ϕm ˙ ϕs) 2 tanh(β ˙ϕm) + Fvϕ˙m. (5.1)

The behaviour of the modelled fault is plotted in Figure 5.3. The figure shows several superimposed friction curves illustrating an increasing fault, where the lowest curve represents the normal friction curve, i.e. the non-faulty case and the uppermost torque curve represents the fault used in this analysis. The model’s parameters Fc, Fs and ϕsare variable and change for every increase in the

mag-nitude of the fault. The variable ˙ϕmdenotes the motor angular velocity. Added

friction torque naturally results in a larger applied motor torque and Figure 5.4 shows plots of the torque and velocity signals when affected by the friction fault and when not affected by the friction fault, on top of each other. The markings in the figure hints where the largest differences appear for the first time in the periodic signal. It should be noted that the model is not derived from a specific robot. It has, nevertheless, authentic characteristics seen from previous robot friction identification situations.

(41)

5.2 Disturbances 29 0 100 200 300 400 0.1 0.15 0.2 0.25

Torque [Nm]

Figure 5.3. Ten torque curves representing increasing friction fault.

0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 Time [s] Torque [Nm]

(a) Effect of friction fault on the torque signal.

0 0.5 1 1.5 −400 −300 −200 −100 0 100 200 300 400 Time [s]

(b) Effect of friction fault on the velocity signal (the fault has no effect on the velocity).

Figure 5.4. Faulty signals superimposed on non-faulty signals.

5.2.2 Shift and delay

Other disturbances that have been represented in the signal model are angle shift and delay. The angle shift is a disturbance which occurs when a robot arm reaches a too high velocity and performs a larger angle displacement than anticipated, resulting in a slight shift in time causing a desynchronisation with nominal mea-surements (of exactly the same test cycle). The delay disturbance occurs when the robot arm moves through a zero-velocity point, also referred to as a fine-point, and remains there for too long and therefore causing a slight shift in time, i.e. a desyn-chronisation with nominal measurements. These disturbances have been seen to

(42)

30 Measurements and disturbance models appear in the practical application of the method on physical robots and are most likely a result of the robot manipulator being a mechanical device and subject to inevitable deviation from one occasion to another (even though performing the same test cycle). The disturbances do not have any analytical representations as the friction torque has but are represented by direct manipulation in the test data. The effects of these disturbances are illustrated in Figure 5.5 and Figure 5.6. The markings in the figures show where the effect of the fault appears for the first time in the periodic signal.

0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 Time [s] Torque [Nm]

(a) Effect of angle shift on torque data.

0 0.5 1 1.5 −400 −300 −200 −100 0 100 200 300 400 Time [s]

(b) The angle shift occurs every time the velocity reaches an extremum point.

Figure 5.5. Faulty signals superimposed on non-faulty signals.

0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 Time [s] Torque [Nm]

(a) Effect of delay on torque data.

0 0.5 1 1.5 −400 −300 −200 −100 0 100 200 300 400 Time [s]

(b) The delay occurs at every zero-velocity point.

(43)

5.3 Residuals 31

5.2.3 Noise

The noise that we expect in the torque and velocity data, which is a simulation of the noise in the resolver, is in this investigation approximated as

nτ(s) = − 4.2s3_{+ 154s}2_{+ 140s} s3_{+ 32s}2_{+ 6.1 · 10}3_{s + 3 · 10}5· u(s) (5.2) nϕ˙(s) = s2 s2_{+ 200s + 10}4 · v(s) (5.3)

for the torque and velocity signal respectively. Here u(s) and v(s) are white noise sequences both with variance ν = (0.5 · 103₎2_{. Figure 5.7 shows the periodograms}

of nτ(t) and nϕ˙(t). 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 Frequency [Hz] Density

(a) Periodogram of the signal nτ(t).

0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 Frequency [Hz] Density

(b) Periodogram of the signal nϕ˙(t).

Figure 5.7. Periodogram of noise models.

5.3 Residuals

The residuals can be considered as generated through the scheme described in Section 3.4.3, i.e. as a subtraction between process and model output, where past (nominal) measurements can be seen as the model output. We will continually work with these two residuals, the torque residual and the velocity residual. Figure 5.8 shows an example of how the respective residuals may look in the time and frequency domain when the signals are affected by the maximum values (used in the analysis) of the disturbances and the friction fault.

(44)

32 Measurements and disturbance models 0 0.5 1 1.5 2 2.5 3 −3 −2 −1 0 1 2 3 Time [s] Torque [Nm]

(a) Torque residual affected by friction fault, de-lay, shift and noise.

0 0.5 1 1.5 2 2.5 3 −1000 −500 0 500 1000 Time [s] Velocity [rad/s]

(b) Velocity residual affected by friction fault, delay, shift and noise.

0 5 10 15 20 0 0.02 0.04 0.06 0.08 0.1 Frequency [Hz] Density

(c) Torque residual in frequency domain.

0 5 10 15 20 0 0.5 1 1.5 2 2.5x 10 4 Frequency [Hz] Density

(d) Velocity residual in frequency domain.

Figure 5.8. Torque and velocity residuals with added noise nτ(t) and nϕ˙(t), in time and frequency domain.

5.3.1 Noise sensitivity

It is always of interest to assess how much of a signal is useful data and how much is noise. This is commonly evaluated by calculating the signal-to-noise ratio (SNR) of the signal. In our case we had the noise models nτ(s) and nϕ˙(s) and

residuals τres(t) and ˙ϕres(t) to calculate the measure. The SNR is for the torque

and velocity residual calculated as

SN Rτ = P τres(k)2 P nτ(k)2 SN Rϕ˙ = P ˙ϕres(k)2 P nϕ˙(k)2 (5.4) where τres(k) and ˙ϕres(k) are the residuals of respective signal. The SNR

com-puted on the residuals shown in Figure 5.8 is

(45)

Chapter 6

Design of experiment

In this chapter the design of experiment used in the analysis will be explained as well as some brief background theory. The design described here was constructed by Dr. Kari Saarinen, principal scientist at ABB in Västerås.

6.1 Introduction

One of the most important steps in research often pointed out by statisticians, is the conduct of efficient and reliable experiments and doing so in a systematic manner. It is in the literature advised not to perform experiments arbitrarily and change control factors one at a time, but to utilise variation in multiple vari-ables simultaneously and to let them do so in a systematic way. This is called constructing factorial experiment designs [1]. Assuming that researchers typically face limited resources in terms of time and money, to use them efficiently, obviously is of benefit.

Terminology

Treatment A controllable factor of influence (on the yield), also re-ferred to as an experiment design variable.

Treatment level A symbol for when a treatment is high, low or in the middle.

Run One experiment trial.

Experiment project A set of runs.

6.2 Fractional factorial designs

When a researcher is about to construct an experiment project (set of runs) he or she needs to decide which influences (treatments) to test the response to. A simple example could be deciding whether to test plant growth with respect to varying sun light or with respect to amount of given water, or both! The central point in factorial design is to assign these factors of influence a discrete set of levels in

(46)

34 Design of experiment Treatments Run a b c 1 - - -2 - - + 3 - + -4 - + + 5 + - -6 + - + 7 + + -8 + + +

Table 6.1. A full factorial matrix.

Treatments Run d a b c 1 - - - -2 + - - + 3 + - + -4 - - + + 5 + + - -6 - + - + 7 - + + -8 + + + +

Table 6.2. Augmented matrix with an

extra factor d by confounding a, b and

c.

which they are let to vary. Levels may for example be represented as (+) and (-) or 1 and 0, for high and low levels. Next task being to assign real values to these levels and be consistent in their application throughout the experiment project. A setup of experimental runs can for example look like Table 6.1, where + and - have been used. The number of test runs is determined by the amount of treatments and treatment levels. For a setup of two-level treatments and three-treatment runs we have for a full factorial design a total of 23= 8 runs.

Assume now that we have the 8-run design shown in Table 6.1. Assume further that we only have time to perform half of these experiment runs and thus would like to fractionalise, i.e. choose a subset of the 8-run design. We would do this, as is usually done, by creating a fourth treatment by confounding the existing ones, e.g. letting d = abc, where the sign of d is a multiplication of the signs of

a,b and c. We would then get a setup shown in Table 6.2. From this setup we

choose the runs where d remains constant (either + or -), which for d = +, gives runs number 2, 3, 5 and 8 as our subset. This choice of fraction is said to have the highest resolution from the fact that the variable d was defined using all the original variables. One could of course also choose a subset in a purely random way. For more on this specific issue and experiment designs [1] is recommended.

6.3 Treatments and treatment levels

In our case we have five factors of influence (treatments) shown in Table 6.3. The chosen variables are selected with the expectation that they would be significant in explaining the variation in the test quantities. Three of the variables have three levels while the other two have two levels. This results in when constructing a full factorial design, a total of 108 experiment runs. The used experiment matrix, however, is fractionalised to 42 runs and is shown in Table 6.5 in the end of this chapter. Table 6.3 and Table 6.4 describe the influence factors, their levels and the levels’ values.

(47)

6.4 Experiment matrix 35 Test cycle variables

Angular velocity ϕ˙ Angle movement ϕ Disturbance factors Friction fault f Delay ∆t Angle shift ∆ϕ

Table 6.3. Factors of influence, treatments.

Level

Variables 1 2 3

x1 f Friction fault - 0 1

x2 ϕ˙ Angular velocity 1.2 rad/s 1.5 rad/s 1.8 rad/s

x3 ϕ Angle movement 10◦ 18◦ 26◦

x4 ∆t Delay time - 0 ms 50 ms

x5 ∆ϕ Angle shift −0.5◦ 0◦ 0.5◦

Table 6.4. Treatment levels and their values.

6.4 Experiment matrix

The construction of the experiment design originated from the variables angular velocity ˙ϕ and angular displacement ϕ, both varying in three levels. Moreover,

there are inherent physical constraints between the velocities that can be reached and the displaced angle. The robot arm may for example not be able to reach a specific (high) velocity before it completes its movement if a too small angle is defined, because of the constraints on the maximal acceleration of the robot arm. Therefore the following was stated; The velocity level ˙ϕ = 3 should correspond

to the highest achievable velocity when performing the angle movement with level

ϕ = 2. This makes the combination [ ˙ϕ, ϕ] = [3, 1] unachievable.

Figure 6.1 shows the experimental region with an X-mark for the unachievable experimental situation. 1 2 3 4 5 6 X 3 2 1 1 2 3 ˙ ϕ ϕ

(48)

36 Design of experiment

For every point in the figure the variables ∆ϕ and ∆t are let to vary between 2 and 3 and 1 and 3 respectively, except for point 6. For the experiment situation 6 all variables remain on level 2. This setup resulted in a design matrix of totally 21 trials, as seen as the first 21 trials in Table 6.5.

The two-leveled friction fault f was added at a later point and led to the doubling of the design matrix to totally 42 trials resulting in the final design matrix. Exp. No f ϕ˙ ϕ ∆t ∆ϕ 1 2 1 1 2 1 2 2 1 1 2 3 3 2 1 1 3 1 4 2 1 1 3 3 5 2 1 3 2 1 6 2 1 3 2 3 7 2 1 3 3 1 8 2 1 3 3 3 9 2 2 1 2 1 10 2 2 1 2 3 11 2 2 1 3 1 12 2 2 1 3 3 13 2 3 2 2 1 14 2 3 2 2 3 15 2 3 2 3 1 16 2 3 2 3 3 17 2 3 3 2 1 18 2 3 3 2 3 19 2 3 3 3 1 20 2 3 3 3 3 21 2 2 2 2 2 22 3 1 1 2 1 23 3 1 1 2 3 24 3 1 1 3 1 25 3 1 1 3 3 26 3 1 3 2 1 27 3 1 3 2 3 28 3 1 3 3 1 29 3 1 3 3 3 30 3 2 1 2 1 31 3 2 1 2 3 32 3 2 1 3 1 Exp. No f ϕ˙ ϕ ∆t ∆ϕ 33 3 2 1 3 3 34 3 3 2 2 1 35 3 3 2 2 3 36 3 3 2 3 1 37 3 3 2 3 3 38 3 3 3 2 1 39 3 3 3 2 3 40 3 3 3 3 1 41 3 3 3 3 3 42 3 2 2 2 2

(49)

Chapter 7

Sensitivity analysis

In this chapter we will perform a sensitivity analysis of the chosen test quantities with respect to different test cycles while assuming that the signals are subjected to a friction fault and two disturbances. This will be done using regression analysis. An analysis of main effects will be performed and presented, followed by a similar analysis of an improved model.

7.1 Test quantities

In our diagnosis method the test quantities that are used as base for decision making are the RMS-measure and model fit-measure mentioned in Section 1.3.3. They are computed on the generated residuals. They are repeated below for convenience. φf it=  1 − q PN k=1(sres(k) − ˆsres(k)) 2 q PN k=1(sres(k) − ¯sres) 2  · 100 (7.1) φRM S = s PN k=1s2res(k) N ,s_P_N k=1s2nom(k) N · 100 (7.2) Where s(k) is either τ (k) or ˙ϕ(k).

The objective now is to run an experiment project (using the experiment ma-trix) for each test quantity and use regression to perform a sensitivity analysis, i.e. to analyse how the quantities vary with respect to the chosen treatments (control factors). We will refer to these projects as

• Project 1 Sensitivity analysis of the model-fit measure when applied on the torque signal.

• Project 2 Sensitivity analysis of the model-fit measure when applied on the velocity signal.

(50)

38 Sensitivity analysis

• Project 3 Sensitivity analysis of the RMS-measure when applied on the torque signal.

• Project 4 Sensitivity analysis of the RMS-measure when applied on the ve-locity signal.

7.2 Calculation of main effects

We will now perform a multivariable regression analysis with the test quantities as response variables. The experiment matrix is regarded as the x-matrix. We state the linear model

Y = β0+ β1x1+ β2x2+ β3x3+ β4x4+ β5x5+ . (7.3)

Refer to Table 6.4 for explanation of the variables corresponding to respective xi.

7.2.1 Results

Table 7.1 - Table 7.4 below show the results from the regression. The first column shows the parameter βicorresponding to experiment variable xi while the second

column shows the least-squares estimate of that parameter. The third column shows the probability of the parameter being insignificant in the model, i.e. the probability of the parameter’s confidence interval encompassing the zero value (recall Section 4.2.2). Project 1 β βˆ P(βi= 0) β0 82.5622 8.8009e-08 β1(f ) 1.2343 0.6940 β2( ˙ϕ) -0.4107 0.8248 β3(ϕ) 1.6430 0.3781 β4(∆t) -0.0714 0.9818 β5(∆ϕ) -3.2437 0.0493

Table 7.1. Regression results of the

model-fit measure when computed on the torque signal.

Project 2 β βˆ P(βi= 0) β0 99.4533 3.0339e-51 β1(f ) -2.5584e-14 1.0000 β2( ˙ϕ) 0.2080 0.0514 β3(ϕ) -0.0697 0.5039 β4(∆t) -0.3244 0.0714 β5(∆ϕ) -0.3795 1.4633e-04

Table 7.2. Regression results of the model-fit measure when computed on the velocity signal.

Table 7.1 shows that the model mean β0lies around 82.56 and that the velocity

variable β2, and the disturbance variables β4 and β5 contribute negatively to the

model-fit measure. The friction fault f however contributes positively, which is reasonable since we expect that a fault would contribute to larger values in the test quantities (recall 1.3.3). It is necessary however to note that all β are showing high P -values even if we haven’t chosen any tolerance level to go by.

Table 7.2 shows on one hand a notably high β0value of 99.45, considering that φf it’s largest value is 100, while on the other hand it shows that in fact almost

(51)

7.2 Calculation of main effects 39 parameter of the friction fault β1 is basically zero, i.e. the variable doesn’t seem

to have any impact on the model-fit measure in this setup. Project 3 β βˆ P(βi= 0) β0 128.6522 6.8211e-05 β1(f ) -1.4747 0.8387 β2( ˙ϕ) -5.3737 0.2147 β3(ϕ) -6.5952 0.1299 β4(∆t) -3.1129 0.6681 β5(∆ϕ) -20.4456 2.7857e-06

Table 7.3. Regression results of the

RMS-measure when computed on the torque signal. Project 4 β βˆ P(βi= 0) β0 80.1256 3.0471e-04 β1(f ) -2.1958e-14 1.0000 β2( ˙ϕ) -2.2209 0.4616 β3(ϕ) -10.3871 0.0013 β4(∆t) 2.6357 0.6050 β5(∆ϕ) -15.1859 1.0198e-06

Table 7.4. Regression results of the RMS-measure when computed on the velocity signal.

Recalling that the RMS-measure gives a measure of the energy ratio between the residual and nominal data, one would expect it to be relatively small, since the residual should be a fraction of the nominal data (even in presence of fault). Nevertheless Table 7.3 shows a significantly large β0, implying a variation of the

quantity around approximately 80. Table 7.4 also shows a large β0 albeit smaller

than in Project 3. Another observation is that in both tables the angle parameter

β3 and angle shift parameter β5are the most significant ones.

7.2.2 Analysis of variance

Table 7.5 - Table 7.8 below describe how well the model performs in each project by displaying the R2_{-statistic and the analysis of variance (ANOVA) for every project.}

The R2_{-statistic and the components of the ANOVA are explained in Section 4.2.1}

and Section 4.2.3. The ANOVA basically shows how much of the sums of squares are attributed to the model parameters and how much that are attributed to error, i.e. SSR and SSE. They are used to estimate the significance of the used model,

which is indicated in the F -value and the probability P(βi = 0). They show

the significance of including the “extra” parameters β1, β2,. . . ,βn compared to

only modelling with the constant β0(which can be seen as the “original” model).

More specifically, the last column shows the probability (P(βi = 0)) of the null

hypothesis H0being true, in this case H0: βi= 0, i = 1, 2, . . . , n. How to compute

this probability is briefly described in Section 4.2.3. Fore more on this matter [3, 4] may be consulted. The second column (DF) displays the degrees of freedom of the model and of the error. The degrees of freedom of the model and error are by convention DFM = n−1 and DFE= N −n−1 respectively, where N is the number

of samples and n is the total number of model parameters. The fourth column (SS/DF) shows the sum of square for respective source averaged over degrees of freedom.

(52)

40 Sensitivity analysis

Looking at Table 7.5 the P -value shows that the used model isn’t that much better than it would be using only the constant β0to model the response variable.

The regression coefficients run a 42% “risk” of being insignificant in modelling the response variable. This is even more evident when looking at the R2-statistic.

Project 1 R2 _0.1239 ANOVA Source DF SS SS/DF F P(βi= 0) Regression 5 517.8910 101.7056 1.0184 0.4213 Error 36 3.6614e+03 103.5782 Total 41

Table 7.5. Residual analysis of Project 1.

We see from Table 7.6 - Table 7.8 that the model on average “catches” about 50% of the variation in the data indicated by the R2_{-value. That’s quite a low}

number, not mentioning its performance in Project 1. Even so it is rather arbitrary to state what’s low and what’s not, there are no specific rules here to follow, and we have yet to decide what can be considered as a sufficient P -value.

Project 2 R2 _0.4153 ANOVA Source DF SS SS/DF F P(βi= 0) Regression 5 8.1720 1.6344 5.1137 0.0012 Error 36 11.5060 0.3196 Total 41

Table 7.6. Residual analysis of Project 2.

Project 3

R2 0.5023

ANOVA

Source DF SS SS/DF F P(βi= 0)

Regression 5 1.9728e+04 3.9457e+03 7.2658 8.5829e-05

Error 36 1.9550e+04 543.0459

Total 41

Test Cycle Optimization using Regression Analysis

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Test Cycle Optimization using Regression Analysis.

Test Cycle Optimization using Regression Analysis.

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Abstract

Acknowledgements

Contents

Chapter 1

Introduction

1.1

Definition of a robot

1.2

Fault Detection in robots

1.3

The diagnosis method

1.3.1

Test cycle

1.3.2

Test data

1.3.3

Test quantities

1.3.4

Diagnosis statement

1.4

Outline

1.5

Limitations

Chapter 2

Robotics

2.1

Introduction

2.2

Homogeneous transformations

2.3

Rigid body modelling

Chapter 3

Fault detection and isolation

3.1

Introduction

3.2

Limit checking

3.3

Hardware redundancy

3.4

Model-based fault detection

3.4.1

Fault modelling

3.4.2

Residuals and residual generation

3.4.3

Fault detection using parity equations

3.4.4

Fault detection using signal models

3.4.5

Observer based fault detection

3.5

Change detection

3.6

Fault isolation

Chapter 4

Modelling of time series

4.1

Introduction

4.2

Linear regression

4.2.1

The R

-statistic

4.2.2

t-statistics

4.2.3

F -statistics

4.3

Autoregressive models

4.3.1