Appendix

Back to Abstract, References.

Matlab Programs to Perform Analysis and Simulation

% Input Constant Generator

% EE599 Thesis Project

% By: _Elliott Reitz

% Date: 10/20/91

% Input Generator Routine, varying parameter a.

 

% Initialization:

period = 10; % period of input waveform sinusoid

Onoise = .0; % Output noise level from the plant

Inoise = .2; % Input noise level from the operator or control

rand('normal'); % Switch to normal zero-mean noise distribution

p1chg = -1.8;

p2chg = 0;

 

% Initialization continued during 1st two steps of process

i=1;

% Servo system to imitate a real plant

% Second order system definition for Y(t)=k/(s^2+as+b)

p1(i)=5;

p2(i)=113;

a(i) = p1(i)+p2(i); % sum of poles

b(i) = p1(i)*p2(i); % product of poles

k(i) = b(i); % defines unity gain

X1(i) = sin(2*pi*t(i)/period)*(1-Inoise)+Inoise*rand(1);

Y(i) = ( k(i)*X1(i) ) ...

/ ( 1/dt^2 + a(i)/dt + b(i) ) + Onoise*rand(1);

X(i,:) = [ X1(i) 0 0 ];

i=2;

p1(i)=5;

p2(i)=113;

a(i) = p1(i)+p2(i); % sum of poles

b(i) = p1(i)*p2(i); % product of poles

k(i) = b(i); % defines unity gain

X1(i) = sin(2*pi*t(i)/period)*(1-Inoise)+Inoise*rand(1);

Y(i) = ( Y(i-1)*(2/dt^2+a(i)/dt) + k(i)*X1(i) ) ...

/ ( 1/dt^2 + a(i)/dt + b(i) ) + Onoise*rand(1);

X(i,:) = [ X1(i) Y(i-1) 0 ];

 

% Loop for the time observed by this process

for i=3:steps,

if fix(i/10)==i/10, disp(i); end

 

% Perform Drift functions

p1(i) = p1(1) + p1chg * t(i)/t(steps);

p2(i) = p2(1) + p2chg * t(i)/t(steps);

a(i) = p1(i)+p2(i); % sum of poles

b(i) = p1(i)*p2(i); % product of poles

k(i) = b(i); % defines unity gain

 

% Input Waveform: Random noise on a sine wave

X1(i) = sin(2*pi*t(i)/period)*(1-Inoise)+Inoise*rand(1);

X1(i) = sin(2*pi*t(i)/period) + Inoise*rand(1);

% Now define the X where X1=X1, X2=Y(i-1), and X3=Y(i-2)

X(i,:) = [ X1(i) Y(i-1) Y(i-2) ];

 

% Using backward difference approach, Y as a difference equation is:

Y(i) = ( Y(i-1)*(2/dt^2+a(i)/dt) - Y(i-2)/dt^2 + k(i)*X1(i) ) ...

/ ( 1/dt^2 + a(i)/dt + b(i) ) + Onoise*rand(1);

 

end

save xy

 

% RLSM Parameter estimation Routine

% EE599 Thesis Project

% By: _Elliott Reitz

% Date: 10/20/91

% RLSM Parameter estimation Routine

 

% Start the Process with steps with no prior system knowledge

i=1;

A = A0;

ESTk(i) = k(i);

ESTa(i) = a(i);

ESTb(i) = b(i);

ESTp1(i) = ( ESTa(i) - sqrt( ESTa(i)^2-4*ESTb(i) ) )/2;

Ym(i) = A(1)*X(i,1) + A(2)*X(i,2) + A(3)*X(i,3);

e(i) = Ym(i)-Y(i);

i=2;

a(i) = a(i-1); % sum of poles

b(i) = b(i-1); % product of poles

k(i) = k(i-1); % defines unity gain

A = A0;

ESTk(i) = k(i);

ESTa(i) = a(i);

ESTb(i) = b(i);

ESTp1(i) = ( ESTa(i) - sqrt( ESTa(i)^2-4*ESTb(i) ) )/2;

Ym(i) = A(1)*X(i,1) + A(2)*X(i,2) + A(3)*X(i,3);

e(i) = Ym(i)-Y(i);

% Loop for the time observed by this process, Tmax

for i=3:steps,

if fix(i/10)==i/10,disp(i);end

if i<N,ii=i;else,ii=N;end

ff=1-1/(ii); % forgetting factor

 

% RLSM Procedure

K = P*X(i,:)'*inv(ff+X(i,:)*P*X(i,:)'); % Kalman Gain

A = A + K * (Y(i)-A'*X(i,:)'); % Estimate

P = ( eye(3) - K * X(i,:) ) * P / ff; % Covariance

% Use estimates to measure the parameters of the servo

ESTk(i) = -A(1) / (A(3) * dt^2);

ESTa(i) = -(2+A(2)/A(3))/dt;

ESTb(i) = -( 1/(A(3)*dt^2) + 1/dt^2 + ESTa(i)/dt );

ESTp1(i) = ( ESTa(i) - sqrt( ESTa(i)^2-4*ESTb(i) ) )/2;

% Use RLSM outputs to calculate the Model and Model error.

Ym(i) = A(1)*X(i,1) + A(2)*X(i,2) + A(3)*X(i,3);

e(i) = Ym(i)-Y(i);

% Calculate statistics on the process to evaluate the model.

mY = mY+(Y(i)-mY)/ii;

mYm = mYm+(Ym(i)-mYm)/ii;

me = me+(e(i)-me)/ii;

S2Y = S2Y+( (Y(i)-mY)^2 - S2Y)/ii;

S2Ym = S2Ym+( (Ym(i)-mYm)^2 - S2Ym)/ii;

S2e = S2e+( (e(i)-me)^2 - S2e)/ii;

F(i) = (S2Y-S2e)/S2e; % Ftab=2.7

 

t2 = i; % - Last step in plot interval (update on end of loop to incl pt.).

end % of loop on i

 

Param = ESTp1;

save estp1 Param

save lsm

 

% Parameter Prediction Routine with an AR(p/skip)

% EE599 Thesis Project

% By: _Elliott Reitz

% Date: 10/20/91

% Parameter Prediction Routine

 

Y = Param; % the output time series to be estimated

 

% Additional Initialization

t=(1:steps)*dt;

N=100; % number of points in estimation

t1=1;

mY=0;

slope=0;

ssi=0;

S2A=0;

S2Y=0;

conf2=0;

p=200; % Length of AR process

skip=3; % distance between points in AR process

n=fix(p/skip); % no. points in AR process

alfa=.1;

P=alfa*eye(n);

I=eye(n);

A=0*(1:n)';

i=p+1;

% initialize constant part of model

X(1) = 1;

% initialize Ym vector

Ym = 0*(1:steps); % give definition to all Ym vector elements

% to provide continuity in graphs.

 

dstart=(steps/2/skip-fix(steps/2/skip))*skip;

swtch=steps/2-dstart;

start=(n*skip+skip); % start where >p points are available

end

shk(1)=1;

 

% Loop for the time observed by this process

for i=start:steps,

if fix(i/skip)==i/skip, % only estimate on i = multiples of skip

disp(i); % show i

 

if (i-start)/skip<N,ii=(i-start+skip)/skip;else,ii=N;end

ff=1-1/(ii+.0001); % forgetting factor

 

% Select estimation or prediction

if i<swtch,

% allow estimation for first swtch samples

for j=1:n-1, X(j+1) = Y(i-j*skip); end

else,

% predict for the rest of the this series by using Ym not Y as input

for j=1:n-1, X(j+1) = Ym(i-j*skip); end

end

 

% RLSM Procedure

K = P*X'*inv(ff+X*P*X'); % Kalman Gain

if i<=swtch, % estimation stage

A = A + K * (Y(i)-A'*X'); % Estimate

end

P = (I - K * X)*P/ff; % Covariance Matrix

 

% Use RLSM outputs to calculate the Model and Model error.

Ym(i) = A'*X'; % Model output % Model/Pred Output

e(i) = Ym(i)-Y(i); % Est/Pred Error

 

% Calculate statistics on the process to evaluate the model.

mY = mY+(Y(i)-mY)/ii;

S2Y = S2Y+( (Y(i)-mY)^2 - S2Y)/ii;

mYm = mYm+(Ym(i)-mYm)/ii;

S2Ym = S2Ym+( (Ym(i)-mYm)^2 - S2Ym)/ii;

me = me+(e(i)-me)/ii;

S2e = S2e+( (e(i)-me)^2 - S2e)/ii;

 

% Calculate 95% confidance interval and predicted error

if i>swtch,

l=(i-swtch)/skip;

if l>n, jmax=n; else, jmax=l; end % length of AR(n)

for j=2:jmax, % calculate shock vector

shk(j)=shk(1);

for k=2:j, shk(j)=shk(j-1)+A(k)*shk(j-k+1); end

end

shksum2 = 0;

for j=1:jmax, shksum2 = shksum2 + shk(j)^2; end % sum of shk(j)^2

S2el = S2del*shksum2;

conf(i) = 1.96*sqrt(S2el); % Confidance intvl

else,

S2del = S2e; % save for pred stg

conf(i) = 1.96*sqrt(S2e); % Confidance intvl

end

 

else, % fix(i/skip)==i/skip = 0 thus hold values

Ym(i)=Ym(i-1);

conf(i)=conf(i-1);

e(i)=e(i-1);

end % of skip-if

 

t2 = i; % - Last step in plot interval (update on end of loop to incl pt.).

end % of loop on i

 

save bo

 

% Parameter Prediction Routine with an ARMA(1,1)

% EE599 Thesis Project

% By: _Elliott Reitz

% Date: 10/20/91

% Parameter Prediction Routine

 

Y = Param; % the output time series to be estimated

 

% Additional Initialization

t=(1:steps)*dt;

N=100; % number of points in estimation

t1=1;

mY=4.6;

mYm=4.6;

S2A=0;

S2Y=0;

conf2=0;

skip=1; % distance between points in AR process

p=2*skip; % Time-span

n=fix(p/skip); % no. points in AR process

alfa=.1;

P=alfa*eye(n);

Pa=alfa*eye(n);

I=eye(n);

A=0*(1:n)';

Th=0*(1:n)';

swtch=steps/2;

% initialize Ym vector

Ym = 0*(1:steps); % give definition to all Ym vector elements

% to provide continuity in graphs.

 

dstart=(steps/2/skip-fix(steps/2/skip))*skip;

swtch=steps/2-dstart;

start=(n*skip+skip); % start where >p points are available

shk(1)=1;

 

% Loop for the time observed by this process

for i=start:steps,

if fix(i/skip)==i/skip, % only estimate on i = multiples of skip

disp(i); % show i

 

if (i-start)/skip<N,ii=(i-start+skip)/skip;else,ii=N;end

ff=1-1/(ii+.0001); % forgetting factor

 

% RLSM Procedure

if i<swtch, % estimation stage

% allow estimation for first swtch samples

X = [ 1 Y(i-skip) ];

K = P*X'*inv(ff+X*P*X'); % Kalman Gain

A = A + K * (Y(i)-A'*X'); % Estimate of AR(1)

At(i) = Y(i)-A(1)-A(2)*Y(i-skip); % Estimate of shock

Ym(i) = mY*(1-A(2))+A(2)*X(2)-Th(2)*At(i-skip)+At(i);

x = [ 1 -t(i) ];

Ka = Pa*x'*inv(ff+x*Pa*x'); % Kalman Gain

Th = Th + Ka * (Y(i)-Th'*x'); % Estimate of MA(1)

elseif i==swtch,

X = [ 1 Ym(i-skip) ];

At(i) = Y(i)-A(1)-A(2)*Y(i-skip); % Estimate of shock

Ym(i) = mYm*(1-A(2))+A(2)*X(2)-Th(2)*At(swtch-skip)+At(swtch);

K = P*X'*inv(ff+X*P*X'); % Kalman Gain

else,

% predict for the rest of the this series by using Ym not Y as input

X = [ 1 Ym(i-skip) ];

Ym(i) = mYm*(1-A(2))+A(2)*X(2)-Th(2)*At(swtch-skip)+At(swtch);

K = P*X'*inv(ff+X*P*X'); % Kalman Gain

l=(i-swtch)/skip;

end

 

P = (I - K * X)*P/ff; % Covariance Matrix

e(i) = Ym(i)-Y(i); % Est/Pred Error

 

% Calculate statistics on the process to evaluate the model.

mY = mY+(Y(i)-mY)/ii;

S2Y = S2Y+( (Y(i)-mY)^2 - S2Y)/ii;

mYm = mYm+(Ym(i)-mYm)/ii;

S2Ym = S2Ym+( (Ym(i)-mYm)^2 - S2Ym)/ii;

me = me+(e(i)-me)/ii;

S2e = S2e+( (e(i)-me)^2 - S2e)/ii;

 

% Calculate 95% confidance interval and predicted error

if i>swtch,

l=(i-swtch)/skip;

for j=2:l, % calculate shock vector

shk(j)=A(2)^(j-1)*(A(2)-Th(2));

end

shksum2 = 0;

for j=1:l, shksum2 = shksum2 + shk(j)^2; end % sum of shk(j)^2

S2el = S2del*shksum2;

conf(i) = 1.96*sqrt(S2el); % Confidance intvl

else,

S2del = S2e; % save for pred stg

conf(i) = 1.96*sqrt(S2e); % Confidance intvl

end

 

else, % fix(i/skip)==i/skip = 0 thus use previous values for calculations

Ym(i)=Ym(i-1);

conf(i)=conf(i-1);

e(i)=e(i-1);

end % of skip-if

 

t2 = i; % - Last step in plot interval (update on end of loop to incl pt.).

end % of loop on i

 

save box

 

% Motor servo comparison of step response

diary motor.log

% Moter servo comparison of step response

% Third order system vs Second order approximation

% By: _Elliott Reitz

% 12/27/91

 

% Moter servo determination of F on poles

% Third order system vs Second order approximation

% By: _Elliott Reitz

% 12/27/91

 

% Servo Constants

Ra = 5.5; % Ohms

La = .025; % H,

Kg = .01; % Gear Ratio

Kol = 5000; % Open loop gain to set p1=5

Kc = 1; % feedback gain

Kt = .25; % V-Sec/Rad or N-m/Amp,

Jm = 1.19E-4; % Kg-m^2,

F = 1.05E-3; % N-m-Sec,

Kw = Kt/La; % N-m/Amp / H,

ta = La/Ra; % H/Ohms,

tm = Jm/F; % Kg-m^2 / N-m-Sec

 

n1=0;

n2=0;

n3=0;

n4=Kt*Kol*Kg/(tm*ta*Ra*F);

NUM = [n1 n2 n3 n4];

d1=1;

d2=(ta+tm)/(ta*tm);

d3=(Ra*F+Kt*Kw)/(tm*ta*Ra*F);

d4=Kc*Kt*Kol*Kg/(tm*ta*Ra*F);

DEN = [d1 d2 d3 d4];

 

[A,B,C,D]=tf2ss(NUM,DEN)

lam=eig(A)

t=0:.01:1;

Y=step(A,B,C,D,1,t);

 

num = [0 560];

den = [1 117 560];

[a,b,c,d]=tf2ss(num,den)

y=step(a,b,c,d,1,t);

 

thp11

diary off

 

% Motor fault simulation routine

% Moter servo determination of F on poles

% Third order system vs Second order approximation

% By: _Elliott Reitz

% 12/27/91

 

% Servo Constants

Ra = 5.5; % Ohms

La = .025; % H,

Kg = .01; % Gear Ratio

Kol = 5000; % Open loop gain to set p1=5

Kc = 1; % feedback gain

Kt = .25; % V-Sec/Rad or N-m/Amp,

Jm = 1.19E-4; % Kg-m^2,

F = 1.05E-3; % N-m-Sec,

Kw = Kt/La; % N-m/Amp / H,

ta = La/Ra; % H/Ohms,

tm = Jm/F; % Kg-m^2 / N-m-Sec

 

n1=0;

n2=0;

n3=0;

d1=1;

 

if exist('Kol0'), % select part of program to run

% Find Kol for dominant pole near 5.0 for tset = .78 sec

i=0;

for j=1:100;

i=i+1;

Koli(i)=100*j;

n4=Kt*Koli(i)*Kg/(tm*ta*Ra*F);

NUM = [n1 n2 n3 n4];

d2=(ta+tm)/(ta*tm);

d3=(Ra*F+Kt*Kw)/(tm*ta*Ra*F);

d4=Kc*Kt*Koli(i)*Kg/(tm*ta*Ra*F);

DEN = [d1 d2 d3 d4];

[A,B,C,D]=tf2ss(NUM,DEN);

lamk(:,i)=eig(A);

p1k(i)=lamk(3,i);

end

end

 

if exist('Kt0'), % select part of program to run

% Find effect of variation in Kt, armiture torque const,

% caused by loss of motor magnet strength.

i=0;

for j = .01:-.0001:.0001, % increment Kti from Kt to .01*Kt

i=i+1;

Kti(i) = Kt*j;

Kw = Kti(i)/La; % Kw is a function of Kt

n4=Kti(i)*Kol*Kg/(tm*ta*Ra*F);

NUM = [n1 n2 n3 n4];

d2=(ta+tm)/(ta*tm);

d3=(Ra*F+Kti(i)*Kw)/(tm*ta*Ra*F);

d4=Kc*Kti(i)*Kol*Kg/(tm*ta*Ra*F);

DEN = [d1 d2 d3 d4];

[A,B,C,D]=tf2ss(NUM,DEN);

lamkt(:,i)=eig(A);

p1kt(i)=lamkt(3,i);

end

end

 

if exist('F0'), % select part of program to run

% Find effect of increased F, friction on motor

% caused by wear on bushings, or dirt in bushings.

i=0;

for j = 1:10:1000 % increment Fi from F to 15000*F

i=i+1;

Fi(i) = F*j;

n4=Kt*Kol*Kg/(tm*ta*Ra*Fi(i));

NUM = [n1 n2 n3 n4];

d2=(ta+tm)/(ta*tm);

d3=(Ra*Fi(i)+Kt*Kw)/(tm*ta*Ra*Fi(i));

d4=Kc*Kt*Kol*Kg/(tm*ta*Ra*Fi(i));

DEN = [d1 d2 d3 d4];

[A,B,C,D]=tf2ss(NUM,DEN);

lamf(:,i)=eig(A);

p1f(i)=lamf(3,i);

end

end

 

if exist('Ra0'), % select part of program to run

% Find effect of increased Ra, armature resistance,

% caused by overheating of motor windings or warn out

% brushes/commutator.

i=0;

for j = 1:10:1000 % increment Fi from F to 15000*F

i=i+1;

Rai(i) = Ra*j;

ta = La/Rai(i); % function of Ra

n4=Kt*Kol*Kg/(tm*ta*Rai(i)*F);

NUM = [n1 n2 n3 n4];

d2=(ta+tm)/(ta*tm);

d3=(Rai(i)*F+Kt*Kw)/(tm*ta*Rai(i)*F);

d4=Kc*Kt*Kol*Kg/(tm*ta*Rai(i)*F);

DEN = [d1 d2 d3 d4];

[A,B,C,D]=tf2ss(NUM,DEN);

lamr(:,i)=eig(A);

p1r(i)=lamr(3,i);

end

end

 

% Plot routines 

% Graph Selection Constants Call graph w/ thp#

% # = 1 % X1, Y, Ym vs t for the servo

% # = 2 % Parameter k, Est k vs t

% # = 3 % Parameter a, Est a vs t

% # = 4 % Parameter b, Est b vs t

% # = 5 % Parameter p1, Est p1 vs t

% # = 6 % Est Parameter and Predicted Parameter vs t

% # = 7 % Kt vs p1

% # = 8 % F vs p1

% # = 10 % Ra vs p1

% # = 9 % Error, conf vs t

 

% Plot of Input and Output and Modelvs t:

plot(t(t1:t2),X1(t1:t2),t(t1:t2),Y(t1:t2),t(t1:t2),Ym(t1:t2)),..

title('Outputs Y, Ym, and Input X1 vs Time'),..

ylabel('X1, Y, Ym'),xlabel('time'),pause;

plot(t(t1:t2),k(t1:t2),t(t1:t2),ESTk(t1:t2)),..

title('Parameter K and Est of K vs time'),..

xlabel('t'),ylabel('K, ESTk'),pause;

plot(t(t1:t2),a(t1:t2),t(t1:t2),ESTa(t1:t2)),..

title('Parameter a and Est of a vs time'),..

xlabel('t'),ylabel('a, ESTa'),pause;

plot(t(t1:t2),b(t1:t2),t(t1:t2),ESTb(t1:t2)),..

title('Parameter b and Est of b vs time'),..

xlabel('t'),ylabel('b, ESTb'),pause; 

plot(t(t1:t2),p1(t1:t2),t(t1:t2),ESTp1(t1:t2)),..

title('Parameter p1 and Estimated p1 vs time'),..

xlabel('t'),ylabel('p1, ESTp1'),pause;

 

% Plot of EST-Param and Predicted Param vs t:

plot(t(t1:t2),Y(t1:t2),t(t1:t2),Ym(t1:t2),t(t1:t2),...

Ym(t1:t2)-conf(t1:t2),t(t1:t2),Ym(t1:t2)+conf(t1:t2)),..

title('Estimated Parameter, Predicted parameter vs Time'),..

ylabel('Par, PrPar, +conf, -conf'),xlabel('time'),pause;

 

% Plot of p1 vs Kt

plot(Kti,p1kt)

title('Motor Constant Kt vs Dominant Pole p1'),..

ylabel('p1'),xlabel('Kt'),pause;

 

% Plot of p1 vs F

plot(Fi,p1f)

title('Motor Friction Constant F vs Dominant Pole p1'),..

ylabel('p1'),xlabel('Friction, F'),pause;

 

% Plot of Error and confidance intervals

plot(t(t1:t2),abs(e(t1:t2)),t(t1:t2),pe(t1:t2),t(t1:t2),conf(t1:t2)),..

title('Model Error, Predicted Error, and Conf Intvl vs t for ESTb'),..

ylabel('e(t), pe(t), conf(t)'),xlabel('time'),pause;

 

% Plot of p1 vs Ra

plot(Rai,p1r)

title('Motor Armature Resistance Ra vs Dominant Pole p1'),..

ylabel('p1'),xlabel('Ra'),pause;

 

% Plot of Input and Output and Model vs t:

plot(t,Y,t,y),..

title('Step response of Servo and Second Order Approximation'),..

ylabel('Y, y'),xlabel('time'),pause;

 

% Abstract - Fault Prediction With Regression Models