Application Software for a Model Based Approach to Tool Wear Estimation by Kourosh Danai and A. Galip Ulsoy Technical Report No. UM-MEAM-86-36 Department of Mechanical Engineering and Applied Mechanics University of Michigan Ann Arbor, MI 48109-2125

;: I

Chapter 1 INTRODUCTION The purpose of this report is to present the application software developed for a model based approach to tool wear estimation. This approach which uses an adaptive observer for tool wear estimation is based on a dynamic state model of tool wear. The application software for this approach is developed to operate on a Digital Equipment Corporation (DEC) LSI-11/23 plus microcomputer. Although, most of the computer code is in FORTRAN IV and can be used on other computers, certain parts of the code are exclusively written for this type of computer. The developed software can be categorized into three: 1. The program to evaluate the model, and compute its parameters for the purpose of tuning the adaptive observer (see Appendix A). 2. The program to implement the adaptive observer. This program consists of both high and low level software. The high level software is in FORTRAN IV and the low level software is in ASSEMBLY language (see Appendix B). 3. The program to plot and print the estimated results (see Appendix C).

The report describing the methodology, is organized in two parts: 1. Part 1 (Chapter 2) describes the strategy used for evaluating the model and computing the necessary data. 2. Part 2 (Chapter 3) describes the adaptive observer.

Chapter 2 MODEL EVALUATION The design of the adaptive observer is based on a dynamic state model of tool wear. Therefore, the suitability of the model for the purpose of on-line tool wear estimation should be determined. We do this by studying the following points: 1. Simulation results using the nonlinear model. 2. Stability of the system as predicted by the model. 3. Controllability of the model as it relates to parameter estimation. 4. Observability of the wear components by cutting force measurement. Once the suitability of the model for tool wear estimation has been verified, the adaptive observer needs to be designed. The design of the adaptive observer requires knowledge about the canonical parameters, and the observed states. Also, for the purpose of the reconstruction of the wear components from the observed states a transformation matrix T has to be approximated in terms of the cutting conditions, the estimated parameters, and time. This approximation is based on the true value of the transforma

tion matrix obtained from the linearization of the nonlinear model at different operating points. The methodology used to compute the required data is discussed here. 2.1 Simulation The simulation of the nonlinear model was performed by a simulation package developed by Leal Lauderbaugh [1]. This package uses a Runge-Kutta 4 algorithm for integration. 2.2 Linearization In order to (i) study the stability, controllability, and observability conditions of the model; (ii) compute the canonical parameters, and states in observer form for the purpose of tuning the adaptive observer; and (iii) compute the transformation matrix as the basis for the regression analysis, the model has to be linearized. The linearization method used is the subject of this section. The nonlinear model = f'(x, u, t) (2.1) y = g'(x, u, t) (2.2) can be linearized about an operating point (x and u), such that -=Ax+bu (2.3) y -- CT X+ D u (2.4) where x represents the vector of state variables and u is the input variable defining the operating point. In order to obtain such a model x and u are divided into n' intervals

such that, Ax-= X (2.5) ni Au = _ (2.6) n' Then, the Taylor series expansion is used to linearize the nonlinear functions fV(x, u, t) and g'(x, u, t) in each interval, such that 1 af' af8C af2 A = n [1 ax+......+ ] (2.7) X=0 X= ZcX X=(n'-1) LX u=O u ZAu u=(n'-1)Au 1 af' af' af' b = +..... (2.8) X=O X=AX X= (n'-1)AX u=O u=Au u=(n — 1) Au 1 8g' 8g' ag' D =, n' [ +l +,......+ ] (2.9) X=o X=ACX X=(n'-1) AX u=O u= Au u=(n'-1)lAu 1 agi gI' 8g' D = +......+ I (2.10) X=O X= AX X=(nI -1) LLX u=O u=AU u=(n'-1)lu where n' (the number of intervals) in Eqns. (2.5) - (2.6) is selected based on the convergence criterion of the method. For this specific nonlinear model, n' was selected to be 200. The partial derivatives in Eqns. (2.7) - (2.10) are given in [2]. Once the linear model is obtained, it can be used to study its stability condition by eigenanalysis. 2.3 Eigenvalues Eigenvalues of a linear model can be obtained from the roots of the model characteristic equation. The characteristic equation can be obtained by using the relationship det[sI - A] = 0 (2.11)

If for our third order linear model (Eqns. (2.3) and (2.4)) the coefficient matrix A is defined as 1ll1 ca12 cs13 A = C21 Ci22 C23 (2.12) aC31 a32 Ca33 Then the characteristic equation of the model will have the form s3 + els2 + e2s + e3 =O (2.13) where el = - (Call + C22 + C33) (2.14) e2 - C11 C122 + a11 C33 + C122 33S - a12 0a21 - a13 a31 - a23 a32 (2.15) and e3 = - cll C22 Ca33 + 11ll Ca23 a!32 + a3!22 0/13 Ca31 + Ca33 Ca12 Ca21C12 Ca23 a31 - a13a 21 32. (2.16) The above characteristic equation (Eq. (2.13)) can then be solved to obtain the eigenvalues of the model.

2.4 Controllability and Observability According to Takahashi [3] the linear model defined by equations (2.3) and (2.4) is controllable when, Det t = Det [b, A b, Ab....., An-lb] O (2.17) and is observable when Det 0 = Det [c,ATc,. AT)2,....... (A 1c] 0 (2.18) The determinant of the controllability matrix, Q4, and the observability matrix, 8, are calculated. 2.5 Discrete-Time Form The continuous-time linear model defined in Eqns. (2.3) and (2.4) can be transformed into discrete-time form such that x(k + 1) = P x(k) + q u(k) (2.19) y(k) = CT x(k) + D u(k) (2.20) where P = exp(A At) (2.21) q = (P- I) A-1 b (2.22) and At is the sampling time.

For computational purposes, P and q in the above equations can be obtained by the power series expansion of the term exp(A At) [3], such that 1 1. P = [I + (A At) + (A At)2 +.... + +.(A At)m'] (2.23) and q = t [I + -(AAt) +!(AAt)2 +.... + (A At)m' ] b (2.24) 2! 3+ (MI + 1)! where m' in the above equations determines the degree of accuracy of the method. Equations (2.21) and (2.22) are used to discretize the continuous-time model. For that, m = 100 is used. 2.6 Canonical Parameters In order to design the adaptive observer, the parameters of the model in observer canonical form should be computed. For that, the model should first be transformed into observer form. The observer form of the model is, xo(k + 1) = P xo(k) + q u(k) (2.25) y(k) = cT xo(k) + D u(k) (2.26) where P and q in the above equations are defined as -al 1 0 bl P -a2 0 1 q-= b2 (2.27) -as 0 0 b3 J

The parameters of the model in canonical form can be obtained from the transfer function of the discrete-time model defined by Eqns. (2.17) and (2.18)). Such a transfer function will have the form Y (z) bo'? —1 bj z — Y(z) =b Ef + D (2.28) U (z) 1 + li z-i D where the parameters ai and bi are the same as in Eq. (2.25). The above transfer function can be obtained from the relationship G(z) - U() cT [I - P]-q + D. (2.29) Now, if P, q, and cT are defined as P11 P12 P13 ql P21 P22 P23, q = 42, andcT = [ 1 C2 c3 ] (2.30) P31 P32 P33 q3 then the ai and bi will have the form a1 = -(P11 + P22 + P33) (2.31) a2 = PllP22 + PllP33 + P22P33 - P12P21 - P13P31 - P23P32 (2.32) a3 = - PllP22P33 + P1P23Ps2 + P22P13P31 + P33P12-P21 - P12 P23 P31 - P13 P21 P32 (2.33)

and bl = clql + C2q2 + C3q3 (2.34) b2 -ql [- C1 (P22 + P33) + C2P21 + C3p3l] + q2 [C P12 - C2 (P11 + P33) + 3 P32] + q3 [C1P13 + C2P23 - C3 (P11 + P22)] (2.35) b3 = q [C1 (P22P33 - P23P32) + C2 (P31P23 - P21P33) + C3(P21P32 - P31P22)] + q2 [C (P32P13 - p12 P33) + C2 (PlP33 - P31 P13) + C3 (P31 P12 - P11 P32)] + q3 [Cl (P12P23 - P13P22) + C2 (P21P13 - P23Pll) + C3 (PllP22 - P21P12)] (2.36) The above equations (Eqns. (2.29) - (2.34)) are used to compute the parameters ai and bi for the purpose of tuning the parameter estimator. 2.7 The Transformation Matrix The transformation matrix T in the relationship x = T xo (2.37) is required for two purposes:

1. To compute xo from x, so that the accuracy of the estimated states by the state observer can be studied. This study is important for the purpose of tuning the state observer. 2. To develop the data base used as the basis for the regression analysis to approximate the transformation matrix in terms of the cutting conditions, the estimated parameters, and time. The transformation matrix T can be obtained through the observability matrices of the models in the original form and the observer form, such that T = O-i e (2.38) where = [c,PPTc, (pT)2c,......., (pT)n-lc ] (2.39) _ [c,, ()......., (T)n-1 ] (2.40) In the above relationships, P and c are as in Eqns. (2.17) and (2.18) and P and c are as in Eqns. (2.23) and (2.24). Equation (2.36) is used to compute the transformation matrix T. Once the transformation matrix is determined, using Eq. (2.35) the states of the model in observer form are determined. 11

Chapter 3 ADAPTIVE OBSERVER The adaptive observer used here consists of the combined application of parameter estimation and state estimation. In this adaptive observer it is assumed that (i) the system parameters vary slowly as compared with the states, and (ii) the initial observation error due to poor initial parameter convergence is tolerable. This adaptive observer has the form [4], -al 1 0. 0 & 91 -&2 0 1 0 g92 o(k + I) = o(k) + u(k) + [y(k) - ](k)] (3.1) - an 0 0 bm gn (k) =[ 1 0.. o ]o(k)+ Du(k) (3.2) where the ai and bi are estimated on-line by the parameter estimator. A recursive least squares parameter estimation algorithm has the general form [4],

-1 )+ P(k - 2)(k- 1)(3.3) &(k) =~(k - 1)+-~ + 1(k -1)TP(k- 2)0(k- 1) 1 P(k - 2)f(k - 1)q(k - 1)TP(k - 2)1 P(k - 1) = [P(k - 2) - + q(k - 1)TP(k - 2)0(k - 1) (3.4) where y(k) is the value of the measured variable y at time t = kA t for k = 0, 1, 2,...... P(kc) is the matrix of estimation gains,, provides exponential data weighting, and P(k) is the parameter estimation error. +(k) is the vector of measured (or known) variables, and O(k) is a vector of parameter estimates. The above algorithm recursively updates the estimated parameter vector 0(k) defined as (k) = [ a,(k) a2(k)..... a(k) (k) (k)..... (k) ] (3.5) for any process whose equations can be written in the form, y(k) = 0(k - 1)T O(k) (3.6) Thus, the process model must be written in a form that is linear in the unknown parameters, which are the elements of the vector 0(k). The above algorithm can be used for different estimation methods (equation error method and output error method), which differ in defining the vector O(k) and parameter estimation error P(k). These methods are discussed in the next sections. 3..1 Equation Error Method In the equation error method, the vector e(k) and the estimation error i>(k) are defined as

(k-1)T-[ -y(k-1) -y(k-2)..... -y(k-n) u(k-1) u(k- 2).... u(k- m) ] (3.7) and (k) - [y(k) - q(k - )TS(k- 1)1 (3.8) 3..2 Output Error Method In the parameter estimation algorithm based on the output error method the vector,(k) and the estimation error P(k) are defined as 0(k-1)T-[ -9(k-1) -y(k-2)..... -y(k-n) u(k-1) u(k- 2).... u(k- m) ] (3.9) where y(k) = c(k- 1)T - (k) (3.10) y(k) = -(k- 1)T 0(k - 1) (3.11) and the estimation error i(k) is defined as P(k) = y(k) - y(k) + [D(ql) - 1] [y(k) - y(k)]. (3.12) where D(q-l) in the above equation is a fixed moving average filter defined as D(q-) =- 1+ dlq-l +.......+ dq-l (3.13)

Appendix A This appendix presents the computer package developed for evaluating the nonlinear model for tool wear estimation, and computing the necessary data for tuning the adaptive observer. The components of this computer package are: DISCRT: Main program which calls the appropriate subroutines. DATA: Subroutine to collect the necessary information about the operating point for linearization, n' (see Eqns. (2.5) and (2.6)), and m' (see Eqns. (2.21) and (2.22)). CALC: Subroutine to compute the linear model through linearization of the nonlinear model. ROOTS: Subroutine to compute the roots of a polynomial such as shown in Eq. (2.13). CNTRL: Subroutine to compute the controllability matrix. OBSERV: Subroutine to compute the observability matrix. DSCCON: Subroutine to discretize the continuous-time model. CANON: Subroutine to compute the parameters of the model in canonical form. TFOIIM: Subroutine to compute the transformation matrix.

SUBPRO: A package of subroutine for matrix and vector manipulation.

Print file "DISCRT" Page 1 C This program C (1) Asks for the operating points at which the nonlinear C model is to be linearized (SUBROUTINE DATA). C (2) Linearizes the model about-the selected operating point C (SUBROUTINE CALC). C (3) Calculates the eigenvalues (SUBROUTINE ROOTS), and C determinants of the controllability (SUBROUTINE CNTRL) C and observability matrices (SUBROUTINE OBSERV) of the C linearized model. C (4) Transforms the model into discrete time form (SUBROUTINE DSCCON). C (5) Calculates the eigenvalues of the discrete-time form (SUBROUTINE ROOTS). C (6) Calculates the observer canonical form of the model (SUBROUTINE CANON). C (7) Calculates the transformation matrix for the transformation C (SUBROUTINE TFORM). C C C IMPLICIT REAL*8 (A-H, O-Z) C REAL*8 LO C REAL*4 Y,S1,S2,S3,S4,S5,S6 C DIMENSION Q(14),P(7),X(3),A(3,3),B(3),C(3),PD(3,3),AC(3), $ QD(3,3),QDB(3),AA(3,3),ACAN(3),BCAN(3),PHAT(3,3),QHAT(3), $ DOBSV(3,3),C1(3),PDC(3),PDT(3,3),CNOBSV(3,3),T(3,3), $ TT(3,3),GG(3,3),ROOTR(3),RQOTI(3),CTRL(3,3),OBSV(3,3), $ AT(3,3),B1(3),AB(3),ALPHA(3),BETA(3) C COMMON/INOUT / IN, I OUT COMMON/COEF/Q, P, LO, V, Fl, D, RAKE, TEMP 1, TEMP2,FORCE C DATA Y/'Y' / DATA IN/5/,IOUT/7/ C C CALL ASSIGN(5,'TT:') C CALL ASSIGN(7,'TT:') C CALL ASSIGN (6,' LP:') C Now collect the necessary data C C 100 CALL DATA(X,N,DELT,MACK) C C C Compute the matrices A,B,C,D C C CALL CALC(N,A,B,C,DD,X,MACK) C C Now compute the eigenvalues of the continuous-time system C CALL ROOTS(A,N,ROOTR,ROOTI) C C Check the controllability of the system C WRITE (IOUT, 201) 201 FORMAT(lX,'Would you like to check the controllability', $' and observability of the model Y/N?!') READ(IN, 300) S1 17

Print file "DISCRT" Page 2 300 FORMAT (A2) IF(S1.NE.Y) GOTO 101 C CALL CNTRL(NA, BB1, AB, CTRL) CALL DETER (N, CTRL, DET) C C Now check the observability of the system C CALL OBSERV(N, A, C, OBSV, Cl, AC, AT) CALL DETER (N, OBSV, DET) C 101 WRITE (IOUT, 230 ) 230 FORMAT(1X,'Would you like to compute the transfer', $' function in the s-plane, Y/N?!') READ(IN, 300) S6 IF (S6.NE.Y) GOTO 102 C CALL CANON(N, A, B, C, ALPHA, BETA, PHAT, QHAT) C C Compute the matrices P and Q for the discrete time system C 102 WRITE (IOUT, 202) 202 FORMAT(2X,'Would you like to transform the model', $' into discrete-time, Y/N?!') READ(IN, 300) S2 IF(S2.NE.Y) GOTO 900 C CALL DSCCON(AB,NPD,QD,GG,AA, QDB,DELT) C C COMPUTE THE ROOTS OF THE DISCRT-TIME SYSTEM C CALL ROOTS (PD,N, ROOTR, ROOTI) C C Now we can compute the a's and b's of the canonical form C WRITE (IOUT, 203) 203 FORMAT (2X,'Would you like to transform the model', $' into canonical form, Y/N?!') READ(IN, 300) S3 IF(S3.NE.Y) GOTO 900 C CALL CANON(NPD, QDB,C, ACAN, BCAN, PHAT, QHAT) C C C Now we can compute the transformation matrix C WRITE (IOUT, 204) 204 FORMAT(2X,'Would you like to compute the transformation', $' matrix, Y/N?!') READ (IN, 300) S4 IF(S4.NE.Y) GOTO 900 CALL TFORM(N, PD,C,DOBSV, C1,PDC, PDT, PHAT,CNOBSV, T, TT,MT,NT, X, XBAR) C 900 WRITE(IOUT, 330) 330 FORMAT(//,2X,'Would you like to try another set of input', $ 2X,' Y/N?!',3X) READ(IN, 300) S5 IF(S5.EQ.Y) GOTO 100 18

Print file "DISCRT" Page 3 C STOP END 19

Print f ile "DATA" Page 1 SUBROUTINE DATA (X, N, DELT, MACK) C C This subprogram collects the input data C IMPLICIT REAL*8 (A-H,O-Z) REAL* 8 LO REAL*4 Y,S1 C DIMENSION Q(14),P(7),X(3) C COMMON/ INOUT /IN, IOUT COMMON/COEF/Q, P LO, V, Fl, D, RAKE, TEMP 1, TEMP2, FORCE C DATA Y/'Y' / C WRITE (IOUT, 100) 100 FORMAT (5X,'Enter the cutting speed, in m/min,',5X) READ (IN, 101) V 101 FORMAT (F10.5) C WRITE (IOUT, 102) 102 FORMAT(5X,'Enter the feed rate, in mm/rev,',5X) READ(IN, 101) F1 WRITE (IOUT, 103) 103 FORMAT (5X,'Enter the depth of cut, in ram,',5X) READ(IN, 101) D C C WRITE(IOUT, 106) C 106 FORMAT(5X,'Enter the rake angle, in degrees',5X) C READ (IN, 101) RAKE C RAKE=RAKE* (2. ODO*3.141593)/360. ODO RAKE=10. ODO* (2. ODO*3.141593)/360. ODO C N=-3 C WRITE (IOUT, 107) 107 FORMAT(5X,'Enter the width of flank wear', $' due to abrasion, in mm', 5X) READ(IN, 101) X(1) C WRITE (IOUT, 109) 109 FORMAT 5X,'Enter the width of flank wear', $' due to diffusion, in mm,',5X) READ(IN, 101) X(2) C WRITE (IOUT, 110) 110 FORMAT(5X,'Enter the depth of crater wear,', $'in mm,') READ(IN, 101) X(3) C WRITE (IOUT, 117) 117 FORMAT(5X,'Enter the number of steps for linearization') READ (IN, 118) MACK 118 FORMAT(I3) C 116 WRITE(IOUT,111) 111 FORMAT(5X,'Enter the time increment for simulation,', $'in min.',3X) READ(IN,112) DELT 112 FORMAT(F8.6) 20

Print file "DATA" Page 2 C WRITE (6,113) V, F1,D,X(1),X(2),X(3),MACK, DELT 113 FORMAT(1OX,'v =',G15.8,/,10X,'f =',G15.8,/, $ 9X,' d =',G15.8,/,7X,' wfl =',G15.8,/, $ 7X,' wf2 -',G15.8,/,8X,' wc -',G15.8,/, $ 6X,' MACK =',I3,/,6X,-' DELT =',G15.8,//) C RETURN END 21

Print file "CALC" Page 1 SUBROUTINE CALC (N, A, B, C, DD, X, MACK) C C C This subprogram computes the coefficient matrices of the linear model C C IMPLICIT REAL*8 (A-H,O-Z) REAL* 8 LO C DIMENSION A(N,N),B(N),C(N),Q(14),P(7),X(N) C COMMON/ INOUT / IN, IOUT COMMON/COEF/Q, P LO, V, F1, D, RAKE, TEMP 1, TEMP2, FORCE C C C C C DO 2 L=1,N DO 2 M=1,N A(L,M) =0.OD0 B(L)=0. OD0 2 CONTINUE C G1=V/LO G2=Q(1) *DCOS (RAKE) /F1 C A(1, 1) =G1* (-1. ODO+G2*Q (13) ) A(1,2)=G1*G2*Q(13) A (1, 3) =-G1*G2*Q (14) /D B(1)=G1*G2* (Q (9) *F1** (P (7) -1.0D0) * (1.D0-Q (10) *RAKE) - $ Q (11)/F1-Q (12) *V/F1) C RACK=DBLE (FLOAT (MACK)) +1. ODO DO 20 I=0,MACK BACK=DBLE (FLOAT (I)) VAR1=BACK*X (1) /RACK VAR2=BACK* X (2) /RACK VAR3=BACK*X (3) /RACK VAR4=BACK*F1 /RACK C C IF(VAR4.GT.0.02D0) GOTO 201 FORCE=0.ODO TEMP 1=0.ODO TEMP2=0.ODO GOTO 202 201 FORCE= (Q(9)*VAR4**P(7) * (1.D0-Q(10) *RAKE) -Q(11) -Q(12) *V) *D+ $ Q(13) *D* (VAR1+VAR2) -Q(14) *VAR3 C TEMP1=Q (6) *V**P (1) *VAR4**P (2) +Q (7) * (VAR1+VAR2) **P (3) C TEMP2=Q (8) *FORCE *V* *P (4) *VAR4**P (5) *D**P (6) C 202 G3=27 3. OD0+TEMP1 G4=273. OD0+TEMP2 G5=DEXP (-Q ( 3 )/G3) G6=Q (2) *DSQRT (V) G7=Q(4) *V G8=DEXP (-Q (5) /G4) IF(VAR4.GT.0.001D0) GOTO 203 22

Print file "CALC" Page 2 G9=0.ODO G10=0. ODO Gll=0. ODO GOTO 204 203 G9=Q(8) *V**P(4) *VAR4**P(5) *D**P (6) G1 —=Q (9) *P (7) *VAR4** (P (7)-1.DO) * (1. ODO-Q (10) *RAKE) *D Gll=Q (8) *FORCE*V**P (4) *P (5) *VAR4** (P (5)-l.DO) *D**P (6) C 204 A (2,1) =A (2, 1) +(G6*Q(3) *Q(7) *P (3) * (VAR1+VAR2) ** (P (3) -1. ODO) *G5/ $ (G3*G3) )/RACK A(2,2)=A(2,1) A(2,3)=0.OD0 IF(VAR4.GT.0.02D0) GOTO 205 B(2)=0.ODO GOTO 206 205 B (2) =B (2) + (G6*G5*Q (3) *Q (6) *V**P (1) *P (2) *VAR4** (P (2) -1.OD0)/ $ (G3*G3) )/RACK C 206 A (3, 1) =A (3, 1) + (G7*G8* (Q (13) *D+ (FORCE*Q (5) *Q (13) *D*G9)/ $ (G4*G4)))/RACK A(3,2)=A(3,1) A(3, 3) =A(3, 3) + (G7*G8* (-Q (14) - (FORCE*Q (5) *Q (14) *G9) / $ (G4*G4)))/RACK B (3)=B (3)+ (G7*G8* (G10+ (FORCE*Q (5)* (G9*G10+Gll))/(G4*G4)))/RACK C C CFLANK=(A(2,1)*VAR1+A(2,2)*VAR2+A(2,3)*VAR3+B(2)*VAR4)/ C $ (BACK+1.ODO) C CCRATE=(A(3,1)*VAR1+A(3,2)*VAR2+A(3,3)*VAR3+B(3)*VAR4)/ C $ (BACK+1.ODO) C 20 CONTINUE C AFLANK=WF1 (X(1),X(2),X(3),F1) DFLANK=WF2 (X(1),X(2),X(3),F1) DCRATE=WC(X(1),X(2),X(3),F1) CAFLAN= (A(1, 1) *X(1)+A(1,2) *X(2)+A(1,3) *X(3)+B(1) *F1) CFLANK= (A(2,1)*X(1)+A(2,2) *X(2)+A(2,3) *X(3) +B (2) *F1) CCRATE=(A(3, 1) *X (1) +A(3,2) *X (2) +A (3, 3) *X (3) +B (3) *F1) WRITE ( 6, 305) AFLANK, CAFLAN, DFLANK, CFLANK, DCRATE, CCRATE 305 FORMAT(' ACTUAL ABRASIVE FLANK RATE=',G16.8,10X, $'COMPUTED ABRASIVE FLANK RATE=',G16.8,/, $' ACTUAL DIFFUSIVE FLANK RATE=',G16.8,10X, $'COMPUTED DIFFUSIVE FLANK RATE=',G16.8,/, $' ACTUAL CRATER RATE=',G16.8,19X, $'COMPUTED CRATER RATE=',G16.8) C C Record the components of the C matrix C C(1)=Q(13) *D C(2)=Q(13)*D C (3)=-Q (14) C DD=(Q (9) *Fl** (P (7)-1. ODO) * (1.DO-Q (10) *RAKE) -Q (11)/F1-Q (12) *V/ $ F1)*D WRITE(6,160) 160 FORMAT(///,28X,'THE A MATRIX',//) DO 15 I=1,N 15 WRITE(6,161) (A(I,J),J=1,N) 161 FORMAT(3(5X,G16.8)) C 23

Print file "CALC" Page 3 WRITE (6, 162) 162 FORMAT(///,28X,'THE B VECTOR',//) DO 16 I=1,N 16 WRITE(6,163) B(I) 163 FORMAT (27X, G1 6.8) C WRITE (6, 164) 164 FORMAT(///, 28X,'THE C VECTOR', //) WRITE(6,165) (C(I),I=1,3) 165 FORMAT(3 (5X, G16.8)) C C WRITE(6,168) 168 FORMAT(///,25X,'THE D COMPONENT',//) WRITE(6,169) DD 169 FORMAT (25X, G16.8) C RETURN END 24

Print file "ROOTS" Page 1 SUBROUTINE ROOTS (A, N,ROOTR, ROOTI) C C PROGRAM FOR REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL C IMPLICIT REAL*8 (A-H, O-Z) C DIMENSION A(N,N),CON(4),W(4),ROOTR(3),ROOTI(3) C COMMON/ INOUT/ IN, IOUT C C C $ A(2,2) *A(1,3)*A(3,1)+A(3,3)*A(1,2)*A(2,1)$ A(1,2) *A(2,3) *A(3, 1) -A (1, 3) *A( 3, 2) C CON(2)= A(1,1) *A(2,2)+A(1,1)*A(3,3)+A(2,2) *A(3,3) - $ A(1,2) *A(2,1)-A(1,3) *A(3,1) -A(2,) *A(3,2) C CON (3) =- (A(1, 1) +A(2,2) +A (3,3)) C CON(4)=1.0D0 C CALL POLRT (CON, W, 3, ROOTR, ROOTI, IER) C IF(IER-1) 90, 60,70 C 60 WRITE (6, 65) 65 FORMAT(////34H ORDER OF POLYNOMIAL LESS THAN ONE) GO TO 100 C 70 IF(IER-3) 75,80,78 C 75 WRITE(6,77) 77 FORMAT(////36H ORDER OF POLYNOMIAL GREATER THAN 36) GOTO 100 C 78 WRITE(6, 79) 79 FORMAT(////31H HIGH ORDER COEFFICIENT IS ZERO) GOTO 100 C 80 WRITE ( 6., 85) 85 FORMAT(////50H UNABLE TO DETERMINE ROOT. THOSE ALREADY FOUND ARE) GOTO 100 C 90 WRITE(6, 95) 95 FORMAT(////5X, 9HREAL ROOT, 6X, 12HCOMPLEX ROOT//) C DO 96 I=1,3 96 WRITE (6, 97) ROOTR (I), ROOTI (I) 97 FORMAT(2E16.7) C 100 RETURN END C SUBROUTINE POLRT C PURPOSE C COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL 25

Print file "ROOTS" Page 2 C C USAGE C CALL POLRT (XCOF,COF,M, ROOTR,ROOTI, IER) C C DESCRIPTION OF PARAMETERS C XCOF -VECTOR OF M+1 COEFFICIENTS OF THE POLYNOMIAL C ORDERED FROM SMALLEST TO LARGEST POWER C COF -WORKING VECTOR OF LENGTH M+1 C M -ORDER OF POLYNOMIAL C ROOTR-RESULTANT VECTOR OF LENGTH M CONTAINING REAL ROOTS C OF THE POLYNOMIAL C ROOTI-RESULTANT VECTOR OF LENGTH M CONTAINING THE C CORRESPONDING IMAGINARY ROOTS OF THE POLYNOMIAL C IER -ERROR CODE WHERE C IER=0 NO ERROR C IER=1 M LESS THAN ONE C IER=2 M GREATER THAN 36 C IER=3 UNABLE TO DETERMINE ROOT WITH 500 INTERATIONS C ON 5 STARTING VALUES C IER=4 HIGH ORDER COEFFICIENT IS ZERO C C REMARKS C LIMITED TO 36TH ORDER POLYNOMIAL OR LESS. -C FLOATING POINT OVERFLOW MAY OCCUR FOR HIGH ORDER C POLYNOMIALS BUT WILL NOT AFFECT THE ACCURACY OF THE RESULTS. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C NEWTON-RAPHSON ITERATIVE TECHNIQUE. THE FINAL ITERATIONS C ON EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL C RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED C ERRORS IN THE REDUCED POLYNOMIAL. C C SUBROUTINE POLRT (XCOF,COF,M, ROOTR, ROOTI, IER) DIMENSION XCOF(1),COF(1),ROOTR(1),ROOTI(1) DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2, SUMSQ, 1 DX,DY, TEMP,ALPHA,DABS C C C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION C STATEMENT WHICH FOLLOWS. C DOUBLE PRECISION XCOF, COF,ROOTR, ROOTI C C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS C ROUTINE. C THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE C CONSTANT IN STATEMENT 78 TO 1.OD-12 AND IN STATEMENT 122 TO C 1.OD-10. THIS WILL PROVIDE HIGHER PRECISION RESULTS AT THE C COST OF EXECUTION TIME IFIT=0 26

Print file "ROOTS" Page 3 N=M IER=0 IF (XCOF (N+1) ) 10,25,10 10 IF(N) 15, 15,32 C C SET ERROR CODE TO 1 C 15 IER=1 20 RETURN C C SET ERROR CODE TO 4 C 25 IER=4 GO TO 20 C C SET ERROR CODE TO 2 C 30 IER=2 GO TO 20 32 IF(N-36) 35,35,30 35 NX=N NXX=N+ 1 N2=1 KJ1 = N+1 DO 40 L=l,KJ1 MT=KJ1-L+1 40 COF (MT) =XCOF(L) C C SET INITIAL VALUES C 45 XO=.00500101 YO=0.01000101 C C ZERO INITIAL VALUE COUNTER C IN=0 50 X=XO C C INCREMENT INITIAL VALUES AND COUNTER C XO=-10. 0*YO YO=-10. 0*X C C SET X AND Y TO CURRENT VALUE C X=XO Y=YO IN=IN+1 GO TO 59 55 IFIT=1 XPR=X YPR=Y C C EVALUATE POLYNOMIAL AND DERIVATIVES 59 ICTO0 60 UX=0.0 UY=0.0 V =0.0 YT=0. 0 XT=1. 0 27

Print file "ROOTS" Page 4 U=COF (N+1 ) IF(U) 65,130,65 65 DO 70 I=1,N L =N-I+1 TEMP=COF (L) XT2=X*XT-Y*YT YT2=X*YT+Y*XT U=U+TEMP *XT2 V=V+TEMP *YT2 FI=I UX=UX+FI *XT*TEMP UY=UY-F I *YT*TEMP XT=XT2 70 YT=YT2 SUMSQ=UX*UX+UY*UY IF(SUMSQ) 75,110,75 75 DX= (V*UY-U*UX) / SUMSQ X=X+DX DY=- (U*UY+V*UX) / SUMSQ Y=Y+DY 78 IF(DABS(DY)+DABS(DX)-1.OD-12) 100,80,80 C C STEP ITERATION COUNTER C 80 ICT=ICT+1 IF(ICT-500) 60,85,85 85 IF(IFIT)100, 90,100 90 IF(IN-5) 50,95, 95 C C SET ERROR CODE TO 3 C 95 IER=3 GO TO 20 100 DO 105 L=1,NXX MT=KJ1-L+1 TEMP=XCOF (MT) XCOF (MT) =COF (L) 105 COF (L) =TEMP ITEMP=N N=NX NX=ITEMP IF(IFIT) 120,55,120 110 IF(IFIT) 115,50,115 115 X=XPR Y=YPR 120 IFIT=0 122 IF (DABS (Y) -1. OD-4*DABS (X)) 135, 125, 125 125 ALPHA=X+X SUMSQ=X*X+Y*Y N=N-2 GO TO 140 130 X=0.0 NX=NX- 1 NXX=NXX- 1 135 Y=0.0 SUMSQ=0.0 O ALPHA=X N=N- 1 140 COF(2)=COF(2) +ALPHA*COF(1) 145 DO 150 L=2,N 150 COF (L+1) =COF (L+1) +ALPHA*COF(L) -SUMSQ*COF (L-l) 28

Print file "ROOTS" Page 5 155 ROOTI (N2) =Y ROOTR (N2) =X N2=N2+1 IF(SUMSQ) 160,165,160 160 Y=-Y SUMSQ=0. 0 GO TO 155 165 IF(N) 20,20,45 END 29

Print file "CNTRL" Page 1 SUBROUTINE CNTRL (N, A, B, B1, AB, CTRL) C C C This subprogram computes the controllability matrix C C IMPLICIT REAL*8 (A-H, O-Z) C DIMENSION A(N,N),B(N),CTRL(N,N),B1(N),AB(N) C C C DO 4 I=1,N B1(I)=0.0 4 B1l(I)=B1l(I)+B(I) C N1=0 C DO 5 I=1,N CALL PLACE(N,N1,B1,CTRL) CALL MATVEC(N,N,A,B1,AB) DO 13 J=1,N B1(J) =0.0 13 B1(J)=B1(J)+AB(J) 5 CONTINUE C WRITE(6,129) 129 FORMAT(///,20X,'THE CONTROLLABILITY MATRIX',//) C DO 6 I=1,N 6 WRITE(6,130) (CTRL(I,J),J=1,N) 130 FORMAT(3(5X,G16.8)) C RETURN END 30

Print file "OBSERV" Page 1 SUBROUTINE OBSERV(N,A,C,OBSV,C1, AC,AT) C C C C This subprogram computes the observabitily matrix C C IMPLICIT REAL*8 (A-H, O-Z) DIMENSION A(N,N),C(N),OBSV(N,N),C1(N),AC(N),AT(N, N) COMMON / INOUT / IN, I OUT DO 7 I=1,N Cl (I)=0.0 7 C1 (I)=C1 (I)+C(I) C CALL TRANS (N N MAT, NAT, A, AT) C N1=0 C DO 8 I=1,N CALL PLACE (N, N1, C1, OBSV) CALL MATVEC (MAT,NAT,AT, C1,AC) DO 8 K=1,N Cl (K) =0.0 8 Cl (K)=C1 (K) +AC (K) C WRITE (6, 131) 131 FORMAT(//,20X,'THE OBSERVABILITY MATRIX', //) C DO 10 I=1,N 10 WRITE(6, 133) (OBSV(I,J),J=1,N) 133 FORMAT(3 (5X, G16. 8)) C CALL SINVAL (OBSV,D) C RETURN END 31

Print file "DSCCON" Page 1 SUBROUTINE DSCCON (A, B, N, PD, QD, GG, AA, QDB, DELT) C C C This subprogram computes the matrices P and Q for the C discrete time linear system C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(NN), B (N),PD(NN), QD (NN), QDB(N), $ AA(N,N),GG(N,N) C COMMON/ INOUT/ IN, IOUT C C Find the largest component of A C C U=DABS (A(1, 1)) C DO 22 I=1,N C DO 22 J=1,N C V=DABS (A(I, J)) C IF(U.LT.V) U=V C 22 CONTINUE C C Now find the value of Qmin C C KMIN=0 C 408 KMIN=KMIN+1 C GMIN=FLOAT (KMIN) C GMA=FLOAT (N) C VAL= (GMA) *U* (DELT) C CHECK=1.DOO/DGAMMA (KMIN) *VAL** (KMIN) *DEXP (VAL) C IF(CHECK.GT.0.001) GOTO 408 KMIN=100 WRITE ( 6, 180) KMIN 180 FORMAT(//,25X,'number of iterations =',1X,I3,//) C DO 23 I=1,N DO 23 J=1,N PD(I,J)=0. 0 QD(I,J)=0.0 GG(I, J)=0. 0 IF(I.NE.J) GOTO 90 PD (I,J)=1.0 QD (I,J)=1.0 GG(I,J)=1.0 90 CONTINUE 23 CONTINUE C C DO 31 M=1,KMIN CALL MATMUL (N,N,N,N,GG, A, AA) DO 24 I=1,N DO 24 J=1,N AA ( I, J) =AA ( I, J) *DELT/M PD (I, J) =PD (I, J) +AA(I, J) QD (I, J) =QD (I,J) +AA(I, J) / (M+1) GG(I,J) =AA(I,J) 24 CONTINUE 31 CONTINUE C 32

Print file "DSCCON" Page 2 C DO 25 I=1,N DO 25 J=1,N 25 QD (I, J) =QD (I, J) *DELT CALL MATVEC(N N, QD, BQDB) C WRITE (6, 101) 101 FORMAT(//,25X,' THE P MATRIX',!//) DO 26 I=1,N 26 WRITE(6, 102) (PD (I, J),J=1, N) 102 FORMAT (3 (G16.8,10X)) C WRITE (6, 105) 105 FORMAT(//,25X,' THE Q MATRIX',//) DO 27 I=1,N 27 WRITE(6,103) QDB(I) 103 FORMAT(25X, G16.8) C C RETURN END 33

Print file "CANON" Page 1 SUBROUTINE CANON(N,PD,QDB,C,ACAN,BCAN, PHAT,QHAT) C C THIS SUBPROGRAM CALCULATES THE PARAMETERS OF THE C CANONICAL FORM C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PD(N,N),QDB(N),ACAN(N),BCAN (N),PHAT(N,N), $ QHAT(N),C(N) C COMMON / INOUT / IN, IOUT C C C FIND THE a's FIRST C C C ACAN(1)=-(PD(1, 1)+PD(2,2)+PD(3,3)) ACAN(2)= PD(1,1)*PD(2,2)+PD(1,1)*PD(3,3)+PD(2,2)*PD(3,3)$ PD(1,2)*PD(2,1)-PD (1,3) *PD(3,1) -PD (2,3)*PD (3,2) C ACAN(3) =-PD (1, 1) *PD (2,2) *PD (3,3) +PD (1, 1) *PD (2,3) *PD (3,2) + $ PD(2,2)*PD(1,3)*PD (3,1)+PD(3,3)*PD(1,2) *PD (2,1)$ PD(1,2)*PD(2,3) *PD(31)PD(1,3)*PD(2,1)*PD (3,2) C C NOW LET'S FIND THE b's C BCAN(1)= C(1)*QDB(1)+C(2)*QDB(2)+C(3)*QDB(3) C BCAN(2) =QDB(1)*(-C(1)*(PD (2,2)+PD(3,3)) +C(2)*PD(2,1)+ $ C(3) *PD(3, 1))+QDB(2)*(C(1) *PD (1,2)-C(2)*(PD(1,1) + $ PD(3,3))+C(3) *PD(3,'2))+QDB(3)*(C(1) *PD (1, 3)+C(2)* $ PD(2,3)-C(3) * (PD (1, 1)+PD(2,2))) C BCAN(3)= QDB(1)*(C(1)*(PD(2,2)*PD(3,3)-PD(2,3)*PD(3,2))+ $ C(2)* (PD (3,1)*PD (2,3)-PD (2,1)*PD (3,3))+ $ C(3)*(PD(2,1)*PD(3,2)-PD(3, 1)*PD(2,2)))+ $ QDB(2)*(C(1)*(PD(3,2)*PD(1,3)-PD(1,2)*PD(3,3))+ $ C(2)* (PD(1, 1)*PD(3,3)-PD(3,1) *PD(1,3))+ $ C(3) * (PD(3, 1)*PD(1,2)-PD(1, 1)*PD(3,2)) )+ $ QDB(3) * (C(1) * (PD(1,2) *PD (2,3)-PD(1,3) *PD(2,2))+ $ C(2) * (PD(2,1)*PD (1,3)-PD(2,3)*PD(1,1))+ $ C(3) * (PD (1, 1) *PD (2,2)-PD (2,1) *PD (1,2))) C C C Now write the a's and b's C WRITE(6,150) (I,ACAN(I),I=1,3) 150 FORMAT(//,3(5X,' ALPHA(',I1,')=',1X,G20.12)) C SUMA=1. ODO+ACAN(1) +ACAN(2) +ACAN(3) WRITE(6,159) SUMA 159 FORMAT(/,2X,'1 + SUM OF ALPHA =',G16.8,/) C WRITE(6,151) (I,BCAN(I),I=1,3) 151 FORMAT(//, 3 ( 6X,' BETA(',I1,')=',1X,G20.12)) C SUMB=BCAN(1) +BCAN(2) +BCAN(3) WRITE(6,160) SUMB 160 FORMAT(/,2X,'SUM OF BETA =',G16.8,/) 34

Print file "CANON" Page 2 C C FORM THE CANONICAL MATRICES PHAT AND QHAT C DO 50 I=1,N DO 51 J=1,N PHAT(I, J) =0.ODO IF(J.EQ.1) PHAT(I,J)=-ACAN(I) L=I+1 IF(J.EQ.L) PHAT(I,J)=1. ODO 51 CONTINUE 50 CONTINUE C DO 55 I=1,N QHAT (I) =BCAN (I) 55 CONTINUE C C WRITE (6,120) C 120 FORMAT(//,35X,' THE "PHAT" MATRIX',//) C C DO 52 I=1,N C 52 WRITE(6,121) (PHAT(I,J),J=1,N) C 121 FORMAT(3(10X, G16.8)) C C WRITE (6,123) C 123 FORMAT(//,35X,' THE QHAT VECTOR',//) C DO 53 I=1,N C 53 WRITE(6,124) QHAT (I) C 124 FORMAT (35X,G16. 8) C C RETURN END 35

Print file "TFORM" Page 1 SUBROUTINE TFORM (N, PD, C, DOBSV, C1, PDC, PDT, PHAT, $ CNOBSV, TTT, MT, NT, X,XBAR) C C C THIS SUBPROGRAM COMPUTES THE TRANSFORMATION MATRIX C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION PD(N,N), C (N), DOBSV(N,N),C(N), PDC (N), PDC (N), $ PDT(N,N), PHAT(NN), CNOBSV(NN), T (NN), TT (NN), $ X(N),XBAR(N) C COMMON / INOUT / IN, IOUT C C NOW COMPUTE W TRANSPOSE C CALL OBSERV(N,PD,C, DOBSV, C1, PDC,PDT) C C INVERSE W TRANSPOSE C CALL MINV (DOBSV, N,DETER, C1, PDC) C C NOW COMPUTE W-HAT TRANSPOSE C C(1)=1.0OD0 DO 70 I=2,N C(I)=0.DO 70 CONTINUE C CALL OBSERV(NPHAT,C, CNOBSV, C1, PDC, PDT) C CALL MATMUL(N, N, N, CNOBSV, DOBSV, TT) C C C TRANSPOSE THE TRANSFORMATION MATRIX C C CALL TRANS (N, N, MT, NT, TT, T) C WRITE (6, 430) 430 FORMAT(//,30X,' THE TRANSFORMATION MATRIX',//) C DO 81 I=1,MT 81 WRITE(6,432) (T(I,J),J=1,NT) 432 FORMAT (3 (8X, G20. 12) ) C C NOW CHECK THE TRANSFORMATION MATRIX C CALL MATMUL (N, N, MT, NT, PD, T, TT) CALL MINV(T,N,DETER,C1, PDC) C WRITE (6, 437) 437 FORMAT(//,30X,'INVERSE OF THE TRANSFORMATION MATRIX',//) DO 83 I=1,MT 83 WRITE(6,432) (T(I,J),J=1,NT) C C CALL MATMUL(MT,NT,N,NT, T,TT, PHAT) C C WRITE(6,435) C 435 FORMAT(//,25X,'THE PHAT MATRIX FROM THE TRANSFORMATION',//) 36

Print file "TFORM" Page 2 C DO 82 I=1,N C 82 WRITE(6,436) (PHAT(I,J),J=1,N) C 436 FORMAT(3 (10X,G16.8)) C C CALL MATVEC (N, N, T, X, XBAR) C WRITE (6, 438) 438 ORMAT(//,25X,'THE TRANSFORMED STATES',//) WRITE(6,439) (XBAR(J),J=1,N) 439 FORMAT (30X, G16.8) C RETURN END 37

Print file "SUBPRO" Page 1 SUBROUTINE MATMUL (MA, NA, MB, NB, A, B, AB) C C C This subprogram multiplies the two matrices and puts the C results in matrix AB C C DOUBLE PRECISION A, BAB C DIMENSION A(MA, NA),B(MB,NB),AB (MA, NB) C DO 1 I=1,MA DO 1 J=1,NB AB (I, J)=0.0D0 DO 1 K=1,NA 1 AB(I, J)=AB(I, J) +A(I, K) *B (K,J) C RETURN END C C SUBROUTINE MATVEC (MA, NA, A, V, W) C C MULTIPLIES THE MATRIX A WITH THE VECTOR V C AND RETURNS THE RESULTS IN VECTOR W C DOUBLE PRECISION A,V, W C DIMENSION A(MA,NA),V(NA),W(MA) C C DO 10 I=1,MA W(I)=O. ODO DO 10 J=1INA W(I)=W(I) +A(I, J) *V(J) 10 CONTINUE RETURN END C C SUBROUTINE TRANS (MA, NA, MAT,NAT, A, AT) C C C This subroutine transforms A and puts the result C in AT C DOUBLE PRECISION A,AT C DIMENSION A(MA,NA), AT (NA,MA) C MAT=NA NAT=MA C DO 2 I=1,MAT DO 2 J=1,NAT AT (I, J) =O.ODO 2 AT(I,J)=AT(I,J) +A(J,I) C RETURN END 38

Print file "SUBPRO" Page 2 C C C C SUBROUTINE MATINV (MA, NA, N, A) C C C C Replaces matrix A by its inverse C DOUBLE PRECISION A,DUM C DIMENSION A(MA,NA) C DO 1 I=1,N DUM=A(I, I) DO 3 J=1,N 3 A(I, J) =A(I, J)/DUM A(I, I) =1. ODO/DUM DO 4 J=1,N IF(I.EQ.J) GOTO 4 DUM=A(J, I) DO 5 K=1,N 5 A(J,K)=A(J,K)-A(I,K) *DUM A(J, I) =-A(I, I) *DUM 4 CONTINUE 1 CONTINUE C RETURN END C C C C C SUBROUTINE PLACE (N, N1, BCHECK) C C DOUBLE PRECISION B,CHECK C DIMENSION B(N), CHECK(NN) C N1=N1+1 DO 3 J=1,N K=N1 CHECK (J, K) =0. ODO 3 CHECK (J, K) =CHECK (J, K) +B (J) C RETURN END C C BLOCK DATA C C C This subprogram initializes the COMMON block C 39

Print file "SUBPRO" Page 3 C IMPLICIT REAL*8 (A-H,O-Z) REAL* 8 LO C DIMENSION Q(14),P(7) C COMMON/COEF/Q, P L0, V, F1, D, RAKE, TEMP, TEMP2,FORCE C DATA Q/5.2D-5,20.OD0,8000.OD0,5.OD0,22000. OD0,72.OD0, $ 2500.0D0,0. 056D0,1960.D0, 0. 57D0, 86.OD, 0.D0,500.OD0, $ 2000.ODO/ C DATA P/0.40D0,0.60D0,1.45D0,0.45D0,-0.55D0,-0. 95D0, 0.76D0/ C DATA L0/500.ODO/ C END C C C C FUNCTION DGAMMA(N) C C THIS FUNCTION CALCULATES THE FACTORIAL OF N C C DOUBLE PRECISION DGAMMA C DGAMMA= 1. DO DO 1 I=1,N DGAMMA=DGAMMA* I 1 CONTINUE C RETURN END C C C C SUBROUTINE DETER (N, A, DET) C C THIS SUBROUTINE COMPUTES THE DETERMINANT OF 3 BY 3 C MATRIX C DOUBLE PRECISION A,DET C DIMENSION A(NN) C $ A(1,2) *(A(2,1)*A(3, 3)-A(2,3) *A(3,1))+ $ A(1,3) * (A(2, 1) *A(3,2) -A(2,2) *A(3,1)) C WRITE (6,105) DET 105 FORMAT(///, 5X,' THE DETERMINANT OF THE ABOVE MATRIX IS' $ 2X,G16.8) RETURN END 40

Print file "SUBPRO" Page 4 FUNCTION WF1 (X1, X2,X3,FEED) C IMPLICIT REAL*8 (A-H,O-Z) REAL* 8 LO C DIMENSION Q(14),P(7) C COMMON/COEF/Q, P, L, V, F1, D, RAKE, TEMP 1, TEMP2, FORCE C FORCED=0. ODO IF(FEED.NE.0.ODO) GOTO 201 GOTO 202 201 FORCED=(Q (9) *FEED**P (7) * ( 1.DO-Q (10) *RAKE) -Q (11) -Q (12) *V) * $ D+Q(13) *D* (X1+X2)- Q(14) *X3 C 202 WF1= (V/L0) * (-X1+Q (1) *DCOS (RAKE) *FORCED/ (FEED*D) ) C C C RETURN END C C FUNCTION WF2 (X1, X2, X3, FEED) C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 LO C DIMENSION Q(14),P(7) C COMMON/COEF/Q, P L, V, F1, D, RAKE, TEMP1, TEMP2, FORCE C FTEMP= Q(6)*V**P(1)*FEED**P(2)+Q(7)*(X1+X2)**P(3) C WF2= Q(2)*DSQRT(V)*DEXP ( -Q(3)/(273.D0+FTEMP)) C C WRITE(6,201) FTEMP,WF2 C 201 FORMAT(' FTEMP=',G16.8,2X,'WF2=',G16.8) RETURN END C C C C FUNCTION WC(X1,X2,X3,FEED) C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 LO C DIMENSION Q(14),P(7) C COMMON/COEF/Q, P, L, V, F1, D, RAKE, TEMP 1, TEMP2, FORCE C IF(FEED.NE.0.OD0) GOTO 100 CTEMP=0. ODO FORCED=0.OD0 GOTO 101 100 FORCED= (Q ( 9 ) *FEED * *P ( 7 ) * (1. D 0 -Q (10) *RAKE) -Q (11) -Q (12) *V) * $ D+Q(13)*D* (X1+X2) - Q (144) *X3 CTEMP=Q(8) *FORCED*V**P (4) *FEED**P (5) *D**P (6) 41

Print file "SUBPRO" Page 5 C 101 WC=Q (4) *FORCED*V*DEXP (-Q (5) / (273.D0+CTEMP) ) C C WRITE(6,201) FORCED, CTEMP, WC C 201 FORMAT(' FORCED=',G16.8,2X,'CTEMP=',G16.8,2X,'WC=',G16.8) C RETURN END C C C C C SUBROUTINE MINV C C PURPOSE C INVERT A MATRIX C C USAGE C CALL MINV(A,N,D,L,M) C C DESCRIPTION OF PARAMETERS C A - INPUT MATRIX, DESTROYED IN COMPUTATION AND REPLACED BY C RESULTANT INVERSE. C N - ORDER OF MATRIX A C D - RESULTANT DETERMINANT C L - WORK VECTOR OF LENGTH N C M - WORK VECTOR OF LENGTH N C C REMARKS C MATRIX A MUST BE A GENERAL MATRIX C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C THE STANDARD GAUSS-JORDAN METHOD IS USED. THE DETERMINANT C IS ALSO CALCULATED. A DETERMINANT OF ZERO INDICATES THAT C THE MATRIX IS SINGULAR. C C SUBROUTINE MINV(A,N,D,L,M) DIMENSION A(1),L(1),M(1) C C C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION C STATEMENT WHICH FOLLOWS. C DOUBLE PRECISION A,D, BIGA, HOLD,DABS C C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS C ROUTINE. C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. ABS IN STATEMENT C 10 MUST BE CHANGED TO DABS. 42

Print file "SUBPRO" Page 6 C C SEARCH FOR LARGEST ELEMENT C D=l. 0 NK=-N DO 80 K=1,N NK=NK+N L (K) =K M(K)=K KK=NK+K BIGA=A (KK) DO 20 J=K,N IZ=N* (J-1) DO 20 I=K, N IJ=IZ+I 10 IF(DABS(BIGA)- DABS(A(IJ))) 15,20",20 15 BIGA=A(IJ) L (K) =I M(K) =J 2 0 CONTINUE C C INTERCHANGE ROWS C J=L(K) IF(J-K) 35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A (KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI) =HOLD C C INTERCHANGE COLUMNS C 35 I=M(K) IF(I-K) 45,45,38 38 JP=N* (I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A (JK) A (JK) =A (JI) 40 A(JI) =HOLD C C DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS C CONTAINED IN BIGA) C 45 IF(BIGA) 48,46,48 46 D=0.0 RETURN 48 DO 55 I=1,N IF(I-K) 50,55,50 50 IK=NK+I A (IK) =A (IK) / (-BIGA) 55 CONTINUE C REDUCE MATRIX DO 65 I=1,N 43

Print file "SUBPRO" Page 7 IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K) 60,65,60 60 IF(J-K) 62,65, 62 62 KJ=IJ-I+K A (IJ) =HOLD*A (KJ) +A (IJ) 65 CONTINUE C C DIVIDE ROW BY PIVOT C KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K) 70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE C C PRODUCT OF PIVOTS C D=D*BIGA C C REPLACE PIVOT BY RECIPROCAL C A(KK) =1. 0/BIGA 80 CONTINUE C C FINAL ROW AND COLUMN INTERCHANGE C K=N 100 K=(K-1) IF(K) 150,150,105 105 I=L(K) IF(I-K) 120,120,108 108 JQ=N* (K-1) JR=N* (I-l) DO 110 J=1,N JK=JQ+J HOLD=A (JK) JI=JR+J A(JK) =-A(JI) 110 A(JI) =HOLD 120 J=M(K) IF(J-K) 100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A (KI) JI=KI-K+J A(KI) =-A (JI) 130 A(JI) =HOLD GO TO 100 150 RETURN END 44

Appendix B This appendix presents the software developed for designing the adaptive observer. This software, which is prepared for different stages of the design process, can be categorized as: 1. The computer package to design the parameter estimation algorithms for a linear model. The components of this package are: LST: Main program which calls different parameter estimation algorithms to operate on a linear model. DPAR1: A double-precision least squares parameter estimation algorithm using a diagonal gain matrix and parameter constraints. DPAR2: A double-precision least squares parameter estimation algorithm using a full gain matrix which can be reset periodically. DPAR3: A double-precision least squares parameter estimation algorithm using a full gain matrix which is updated with a constant trace. DPAR4: A double-precision least squares parameter estimation algorithm with U-D factorization of the gain matrix. 45

2. The computer package to apply the adaptive observer to the nonlinear model. For the simulation of the nonlinear model the simulation package SIMULA [1] is used. The components of this package are: LSTNON: Main program to define the nonlinear model and call the necessary subroutines containing the adaptive observer. PARAMI1: A single-precision least squares parameter estimation routine using a diagonal gain matrix and parameter constraints. STATE3: A linear third order state observer. BACK: Subroutine to approximate the transformation matrix from the equations obtained from regression analysis (see [2]), and transform the observed states into estimated wear components. 3. The computer package to implement parameter and state estimation techniques on the measured variable. The components of this package are: SMAIN: Main program to estimate the parameters of the system on-line. SAMPLE: Main program to sample the measured data and store in a data file. ESTIM: Main program to test the tool wear estimation method on the recorded data. SINPUT: Subroutine to collect the necessary data for sampling, such as: (i) sampling period; (ii) the input value; (iii) the voltage range on the analog to digital converter (ADC); (iv) the channel of the ADC; (v) the sampling gain. OPAR21: An output error estimation routine for a first order model. 46

KDADC: Main program to test the ADC board. ADC: Subroutine to operate the ADC (in ASSembly). SSUBS: A package of subroutines to compute the ADC gain, compute the real value of the sampled variable, and reset the real time clock. There are also routines to operate the reat time clock and the digital I/O port which are discussed in detail in [5]. 47

Print file "LST" Page 1 C This program checks the precision of the estimation algorithm C C by C Kourosh Danai C IMPLICIT REAL*8 (A-H,O-Z) REAL* 4 BIT C DIMENSION THETA(6),PHI (6),FGAIN(6,6) DIMENSION BIT(6) C REAL* 4 YESANSWER DATA YES/' Y' / C C INTEGER HOUR,MINUTE, SECOND,MILSEC C C ********** DELETE THIS SET OF CONTROL EQUNS. AND ENTER YOUR'S **** C C NOW SET UP THE MATRICES C C ********** DETERMINE THE VALUE OF N, DCVAL ********************** C N=6 DCVAL=0.2D0 C C ******************************************************************* C CALL ASSIGN(1,' KDLIN.DAT' ) CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') CALL ASSIGN(7,'LP:') C C MINUTE=O C SECOND=O C MILSEC=0O C CALL STRCLK C C C SET THE THETA AND GAIN MATRICES C C 500 M=O WRITE(6,121) 121 FORMAT(' ENTER S1',/) READ(5,122) S1 122 FORMAT(G10.5) C WRITE(6,123) 123 FORMAT(' ENTER S2',/) READ(5,122) S2 C Y=0.ODO YPAST1=0.ODO YPAST2=0. ODO YPAST3=O.ODO UPAST1=0. ODO UPAST2=0. ODO UPAST3=0. ODO 11=0 I2=0 48

Print file "LST" Page 2 C DO 10 I=1,N 10 PHI(I)=O.ODO DO 11 I=1,N DO 11 J=1,N FGAIN(I, J) =O.ODO FGAIN (I, I) =1000.ODO C IF(I.LE.3) FGAIN(I,I)=100.0D0 11 CONTINUE C THETA (1) =-2. 80D0 THETA(2)= 2.80DO THETA (3) =-1. ODO THETA(4)= 8.0ODO THETA(5) =-18. ODO THETA(6)= 7.0DO C A1 —2. 98292630780DO A2= 2.96584835305D0 A3=-0. 982922045614D0 C B1= 10.5360079769DO B2=-21. 0711820226DO B3= 10.5351746106D0 C WRITE(6,125) A1l,A2,A3, B1, B2, B3 D WRITE (7,125) A1,A2,A3,B1,B2,B3 125 FORMAT(/,' TRUE PARAMATERS=', 3(2X, G16. 8 ),/,17X, 3 (2X, G16. 8),///) C C LET' S CALCULATE E(K+1) C 230 Y= A1*YPAST1+A2*YPAST2+A3*YPAST3+ $ B1*UPAST1+B2*UPAST2+B3*UPAST3 C C C YPAST3=YPAST2 YPAST2=YPAST1 YPAST1=-Y C CALL DPAR4 (M,S1, S2, Y, U, FGAIN, PHI,THETA, ERROR) C C LET'S CALCULATE THE NEW INPUT C IF(M.LT.5) GOTO 200 U=DCVAL K=M/ 5 I=K/2 J=I*2 IF(J.EQ.K) U=DCVAL*(1.50DO) L=I/2 J=L*4 IF(J.EQ.K) U=DCVAL* (0.50DO) GOTO 350 200 U-DCVAL 350 UPAST3=UPAST2 UPAS2UPAST2=UPAST1 UPAST1=U C NOW UPDATE THE PHI MATRIX 49

Print file "LST" Page 3 C YTEST=Y/S 1 UTEST=U* S2 CALL SETUP (N, PHI,PHI, YTEST,UTEST) C C STOP THE CLOCK C C CALL STPCLK C WRITE (IOUT, 100) MINUTE, SECOND,MILSEC C100 FORMAT(' TIME IS:',I2,'.',I2,'',I5) C C PREDICT THE OUTPUT C DO 28 I=1,6 BIT(I)=SNGL(THETA(I )) 28 CONTINUE WRITE (1) BIT (1),BIT(2),BIT(3),BIT (4),BIT(5),BIT(6) C IF(M.LE.200) GOTO 230 C WRITE (6, 128) 128 FORMAT(' Would you like to rerun the program?!, Y/N?',/) READ (5, 129) ANSWER 129 FORMAT (A2) IF(ANSWER.EQ.YES) GOTO 500 C CLOSE (UNIT=1) C STOP END 50

Print file "DPAR1" Page 1 SUBROUTINE DPAR1 (MSTEP,S1,S2,Y,U,FGAIN,PHI,THETA,ERROR) C C This is a double precision least squares parameter estimation routine C using a diagonal gain matrix and parameter constraints C IMPLICIT REAL*8 (A-H,O-Z) C ******************************************************* COMMON/ INOUT / IN, IOUT C DIMENSION THETA(6),PHI(6),FGAIN(6,6) DIMENSION BETA(6),C(6),V(6) DIMENSION PAT1 ( 6),PAT2 ( 6), PAT3(6) C C ****************************************************** REAL* 8 LAMDA1, LAMDA2 C ****************************************************** DATA LAMDA1/ 0.95D0/, LAMDA2 /1.0D0 / C N=6 C DO 10 I=1,3 PAT3 (I) =THETA(I) *S1 10 CONTINUE DO 11 I=4,6 PAT3 (I)=THETA(I) /S2 11 CONTINUE C C LET'S CALCULATE E(K+1) C MSTEP=MSTEP+1 CALL VVSCAL(N, PAT3,PHI, Z1) ERROR=Y- Z 1 C C UPDATE D,PAT3 C DO 30 J=1,N PAT1 (J) =FGAIN (J, J) *PHI (J) 30 CONTINUE C DO 31 J=1,N K=J-1 IF(K.EQ.0) GOTO 500 BETA (J) =BETA (K) +PHI (J) *PAT1 (J) GOTO 501 500 BETA (J) =LAMDA1+PHI (J) *PAT1 (J) 501 V(J)=PAT1(J) 31 CONTINUE C C Update the gain matrix C DO 32 J=1,N FGAIN (J, J) =FGAIN (J, J) -PAT1 (J) *PAT1 (J) /BETA (N) 32 CONTINUE C DO 33 I=1,N PAT2 (I) =V(I)/BETA(N) 33 CONT INUE C C UPDATE THETA C 51

Print file "DPAR1" Page 2 DO 13 I=1,N PAT3 (I) =PAT3 (I) +PAT2 (I) *ERROR 13 CONTINUE C C Impose the constraints C GOTO 352 DSET1=-3.10DO*S1 USET1=-2.90DO*S1 IF(PAT3 (1).LT.DSET1.OR.PAT3 (1).GT.USET1) PAT3(1)=-3.0D0*S1 DSET2=2.90D0*S1 USET2=3.10DO*S1 IF(PAT3(2).LT.DSET2.OR.PAT3(2).GT.USET2) PAT3(2)=3. ODO*S1 DSET3=-1.20D0*S1 USET3=-0.80DO*Sl IF(PAT3(3).LT.DSET3.OR.PAT3(3).GT.USET3) PAT3(3)=-1.0D O*S1 USET5=-1.5DO*PAT3 (4) DSET5=-2.5DO*PAT3 (4) IF(PAT3(5).LT.DSET5.OR.PAT3 (5).GT.USET5) PAT3(5)=-2. ODO *PAT3(4) USET6=1.5DO *PAT3 (4) DSET6=0. 5DO *PAT3 (4) IF(PAT3(6).LT.DSET6.OR.PAT3(6).GT.USET6) PAT3(6)=PAT3(4) C.C RESCALE THE THETA VECTOR C 352 DO 21 J=1,3 21 THETA (J) =PAT3 (J) /S1 DO 22 I=4, 6 22 THETA (I) =PAT3 (I) *S2 C C PRINT OUT THE VALUES C WRITE ( 6, 101) MSTEP, Y, ERROR, U D WRITE(7,101) MSTEP,Y,ERROR,U 101 FORMAT(,///,' MSTEP=', I4,/,' Y=' G16.8,10X,'ERROR=',G16.8, $' U=',G16.8) C WRITE ( 6, 102) (THETA(K),K=1, 6) D WRITE(7,102) (THETA(K),K=1, 6) 102 FORMAT(' THETA=', 3(lX,G16.8),/,7X,3(1X,G16.8),//) C RETURN END 52

Print file "DPAR2" Page 1 SUBROUTINE DPAR2 (MSTEP, S 1, S2, Y U FGAIN, PHI, THETA, ERROR) C C This subroutine is a LST with covariance resetting and forgetting factor C estimation algorithm C C by C Kourosh Danai C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION THETA(6),PHI(6),FGAIN(6, 6) DIMENSION PAT1(6),PAT2(6),PAT3(6),PAT4(6),PAT5(6,6) C REAL* 8 LAMDA1, LAMDA2 DATA LAMDA1/0. 95D0 /, LAMDA2/1. ODO / C N=6 C DO 10 I=1,3 PAT3 (I)=THETA(I) *S1 10 CONTINUE DO 11 I=4,6 PAT3 (I) =THETA(I)/S2 11 CONTINUE C C LET'S CALCULATE E(K+1) C MSTEP=MSTEP + 1 C CALL VVSCAL(NPAT3,PHI, Z1) CALL VECMAT (NNPHI, FGAIN, PAT1) CALL VVSCAL(N, PAT1,PHI, Z2) ERROR= (Y-Z1) / (1.0+Z2) C C PRINT OUT THE STUFF C C C C UPDATE THETA C CALL MATVEC (N,N, FGAIN, PHI, PAT2) DO 13 I=1,N 13 PAT3 (I) =PAT3 (I) +PAT2 (I) *ERROR C C UPDATE FGAIN C CALL VVMAT (NNPAT2,PAT1,PAT5) DO 14 I=1,N DO 14 J=1,N FGAIN (I, J)= (FGAIN (I, J)-PAT5 (I,J) / (LAMDA1/LAMDA2+Z2))/LAMDA1 14 CONTINUE C C C RESET THE GAIN MATRIX C GOTO 520 KGAIN=M/10 53

Print file "DPAR2" Page 2 IGAIN=KGAIN/2 JGAIN=IGAIN*2 IF (JGAIN.NE.KGAIN) GOTO 520 DO 24 L=1,N DO 24 J=1,N FGAIN(L, J) -0.ODO0 FGAIN (L, L) =1000. ODO 24 CONTINUE C C NOW BACK TRANSFER THE VALUES C 520 DO 21 J=1,3 THETA(J) =PAT3 (J) /S1 21 CONTINUE DO 22 K=4,6 THETA (K) =PAT3 (K) *S2 22 CONTINUE C WRITE(6,102) MSTEP,Y,ERROR, (PHI(I), I=, 6) D WRITE(7,102) MSTEP,Y,ERROR, (PHI(I),I=1,6) 102 FORMAT(' M=',I4, /,' Y=',G16.8,10X,'ERROR=',G16.8,/, $' PHI =',3(2X,G16.8),/,7X,3(2X,G16.8),/) WRITE(6,104) (THETA(K), K=1,6) D WRITE(7,104) (THETA(K),K=1,6) 104 FORMAT(2X,' THETA =', 3(lX,G16. 8),/,9X, 3(lX, G16.8),//) C RETURN END 54

Print file "DPAR3" Page 1 SUBROUTINE DPAR3 (MSTEP,S 1, S2,Y,U, FGAIN, PHI, THETA, ERROR) C C This subroutine is a double precision constant trace LST C parameter estimation algorithm C C by C Kourosh Danai C C ************************************************************* IMPLICIT REAL*8 (A-H, O-Z) C ************************************************************* C DIMENSION THETA(6),PHI(6),FGAIN(6,6) DIMENSION PAT1(6),PAT2(6),PAT3(6),PAT4(6),PAT5(6,6) C C ************************************************************* REAL*8 LAMDA1 C ************************************************************* DATA LAMDA1/1.ODO/,ALPHA/0.8DO/,DTRACE/10000.ODO/ C N=6 C DO 0. I=1,3 PAT3(I)=THETA(I) *S1 10 CONTINUE DO 11 I=4,6 PAT3(I)=THETA(I) /S2 11 CONTINUE C C LET'S CALCULATE E(K+1) C MSTEP=MSTEP+1 C CALL VVSCAL(N,PAT3,PHI, Z1) CALL VECMAT(N,N,PHI,FGAIN,PAT1) CALL VVSCAL(N,PAT1,PHI, Z2) ERROR=(Y-Z1) / (1.ODO+Z2) C C UPDATE THETA C CALL MATVEC(N,N,FGAIN,PHI,PAT2) DO 13 I=1,N 13 PAT3(I)=PAT3(I) +PAT2 (I) *ERROR C C UPDATE FGAIN C CALL VVMAT(N,N,PAT2,PAT1,PAT5) DO 14 I=1,N DO 14 J=1,N FGAIN(I,J)=FGAIN(I, J) -PAT5 (I,J) / (ALPHA+Z2) 14 CONTINUE C C C RESET THE GAIN MATRIX C GOTO 520 KGAIN=M/10 IGAIN=KGAIN/2 JGAIN=IGAIN*2 IF(JGAIN.NE.KGAIN) GOTO 520 DO 24 L=1,N 55

Print file "DPAR3" Page 2 DO 24 J=1,N FGAIN(L, J)=0.0ODO FGAIN(L,L) =1. ODO 24 CONTINUE C 520 TRACE=O. ODO DO 18 I=1,N TRACE=TRACE+FGAIN ( I, I ) 18 CONTINUE C LAMDA1=TRACE /DTRACE IF(LAMDA1.LE.. ODO) LAMDA1=1. ODO/DTRACE IF (LAMDA1.GT.1.ODO) LAMDA1=1. ODO DO 19 I=1,N DO 19 J=1,N FGAIN (I, J) =FGAIN (I, J) /LAMDA1 19 CONTINUE C C NOW BACK TRANSFER THE VALUES C DO 21 J=1,3 THETA (J) =PAT3 (J)/S1 21 CONTINUE DO 22 K=4,6 THETA (K) =PAT3 (K) *S2 22 CONTINUE C WRITE(6,102) MSTEP,Y, ERROR, (PHI(I),I=1, 6) D WRITE(7,102) MSTEP,Y, ERROR, (PHI (I),I=1, 6) 102 FORMAT(' M=', I4, /,' Y=',G16.8,10X,'ERROR=',G16.8,/, $' PHI =', 3(2X,G16.8),/,7X, 3(2X,G16.8),/) C WRITE(6,103) (THETA(K),K=1, 6) D WRITE(7,103) (THETA(K),K=1, 6) 103 FORMAT(2X,'THETA =',3(lX,G16. 8),/,9X, 3(lX,G16.8),//) C RETURN END 56

Print file "DPAR4" Page 1 SUBROUTINE DPAR4 (MSTEP, S1, S2, Y, U, FGAIN, PHI, THETA, ERROR) C C This subprogram is a double precision LST routine with C U-D decomposition of the gain matrix C C by KOUROSH DANAI C IMPLICIT REAL*8 (A-H,O-Z) C *************************************************************** C DIMENSION THETA(6),PHI(6),FGAIN(6,6) DIMENSION WNEW(6, 6),WPAST (6, 6) DIMENSION BETA(6),WT(6, 6),C(6),V(6) DIMENSION PAT1(6),PAT2 (6),PAT3 (6),PAT4 (6) C REAL* 8 LAMDA1, LAMDA2 DATA LAMDA1/ O. 9 5D0 /, LAMDA2 /1.0D0 / C N=6 C DO 10 I=1,3 PAT4 (I) =THETA (I) *S1 10 CONTINUE DO 11 I=4,6 PAT 4 (I) =THETA (I) / S2 11 CONTINUE C IF(MSTEP.GT.0) GOTO 240 DO 12 J=1,N DO 12 K=1,N IF(J-K) 201,202,203 201 WPAST (J, K) =0.ODO WNEW (J, K) =0.ODO GOTO 12 202 WPAST (J, K) =1 O0D0 WNEW (J, K) =1. ODO GOTO 12 203 WPAST (J, K) =0. ODO WNEW(J, K) =O. ODO 12 CONTINUE C C LET'S CALCULATE E(K+1) C 240 MSTEP=MSTEP+1 CALL VVSCAL (N, PAT4,PHI, Z1) ERROR=Y-Z 1 C C UPDATE UD,PAT3 C CALL TRANS (N,N,WPAST,WT) CALL MATVEC (N, N, WT, PHI, PAT1 ) CALL MATVEC (N, N, FGAIN, PAT1, PAT2) DO 30 J=1,N K=J-1 IF(K.EQ.0O) GOTO 500 BETA (J) =BETA (K) +PAT1 (J) *PAT2 (J) FGAIN(J, J) =BETA(K) *FGAIN(J,J) / (BETA(J) *LAMDA1) C (J) =-PAT1 (J) /BETA (K) 57

Print file "DPAR4" Page 2 GOTO 501 500 BETA(J) =LAMDA1+PAT1 (J) *PAT2 (J) FGAIN (J, J)=LAMDA1*FGAIN (J, J) / (BETA(J) *LAMDA1) C (J) =-PAT1 (J) /LAMDA1 501 V(J)=PAT2(J) C IF(K.EQ.0) GOTO 30 C DO 31 I=1,K WNEW(I, J) =WPAST (I, J) +V(I) *C (J) V(I)=V(I) +WPAST (I, J) *V(J) 31 CONTINUE 30 CONTINUE C C Update the U matrix C DO 32 I=1,N PAT3 (I) =V(I) /BETA (N) DO 32 J=1,N WPAST (I, J) =WNEW (I, J) 32 CONTINUE C C UPDATE THETA C DO 13 I=1,N 13 PAT4 (I)=PAT4 (I) +PAT3 (I) *ERROR C C RESCALE THE THETA VECTOR C DO 21 J=1, 3 21 THETA (J) =PAT4 (J) /S1 DO 22 I=4,6 22 THETA (I) =PAT4 (I) *S2 C C PRINT OUT THE VALUES C WRITE ( 6, 101) MSTEP, Y, ERROR, U D WRITE(7,101) MSTEP, Y, ERROR, U 101 FORMAT(,///,' MSTEP=', I4, /,' Y=', G16.8,10X,' ERROR=', G16.8, $' U=',G16.8) C C WRITE (6, 110) (PHI (J),J=1, 6) C WRITE(7,110) (PHI(J),J=1,6) C 110 FORMAT(' PHI=' 3(2X,G16.8), /,5X, 3(2X,G16.8),/) C WRITE(6, 102) (THETA(K),K=1, 6) D WRITE(7,102) (THETA(K),K=1l,6) 102 FORMAT(' THETA=',3(1X, G16.8),/,7X, 3(1X, G16.8),//) C RETURN END 58

Print file "LSTNON" Page 1 C General Description: This is an adaptive observer routine C for the nonlinear model C C External References: SIMULA C C by KOUROSH DANAI C C This program simulates the tool- wear nonlinear model C C Declarations C LOGICAL*1 ANSW,NO, YES C C ********** CHANGE DIMENSIONS HERE ********* C C THE DIMENSION OF X AND XP =NEQN C THE DIMENSION OF IXC =NCEQN C THE DIMENSION OF Y,IYC AND IREF =NOUT C THE DIMENSION OF U AND IUC =NIN C THE DIMENSION OF IUC =NOUT C THE DIMENSION OF WORK = 5*NEQN C REAL*4 U(1) DIMENSION X(8),XP(8),Y(1),IXC(9),IYC(1),IUC(1),IREF(1) DIMENSION WORK(40) C C Common blocks C COMMON / INOUT / IN, I OUT COMMON/ADCDAC/ADCGAN, IADCOF, IMXADC, DACGAN, IDACOF, DACMIN, DACMAX COMMON/S IZE/NEQN, NIN, NOUT, NCEQN COMMON/CLOCK/H, TEND,TPRINT,DELT,DELDST,DELAY COMMON/INDEX/ICNT, IX, IY, IU, INXC, INYC, INUC, INDV COMMON/MXMN/DVMAX, DVMIN,RIDVMX,RIDVMN C C Data statments C CALL ASSIGN(5,'TT:' ) CALL ASSIGN(6,'TT:' ) CALL ASSIGN(7,'LP:' ) DATA IN/5/,IOUT/6/ C C C set up constants for simulation C NO=' N' YES=' Y' C C ********** CHANGE SIZES HERE ********** C C NIN=# OF PLANT INPUTS C NIN-1 C C NOUT=# OF PLANT OUTPUTS C NOUT=1 C C NEQUN= # OF PLANT STATE EQUATIONS C NEQN=8 59

Print file "LSTNON" Page 2 C C NCEQN= # OF CONTROLLER STATE EQUATIONS C NCEQN=9 C C ********** EDITING IN MAIN COMPLETE ********* C 10 CALL SIMULA(X,XP,Y,U, IXC, IYC, IUC, IREF,WORK) C C Check for a repeat C 9990 WRITE(IOUT,9991) 9991 FORMAT(' Do you wish to repeat the calculations? (y/n)') READ (IN, 9992)ANSW 9992 FORMAT(Al) IF(ANSW.EQ.YES)GO TO 10 IF(ANSW.NE.NO)GO TO 9990 STOP END C SUBROUTINE F(T,X,XP,U) C ********** DO NOT CHANGE THESE DIMENSIONS ********* REAL*4 T,X(1),XP(1),U(1) C ********** DELETE THIS SET OF PLANT EQUNS. AND ENTER YOUR'S **** C COMMON / INOUT / IN, IOUT C C ****** PUT IN THE VALUE OF THE INPUTS ********** C IF(T.GE.0.01) GOTO 200 C AO=500.0 A1=5.2E-05 A2=20.0 A3=8000.0 A4=5.0 A5=22000.0 A6=1960.0 A7=0.57 A8=86.0 A9=0.1 A10=500.0 A11=2000.0 A12=72.0 A13=2500.0 A14=0. 056 C Z1= 0.76 Z2= 0.4 Z3= 0.6 Z4= 1.45 Z5= 0.45 Z6=-0.55 Z7=-0.95 GAM=0. 1745 V-=250. D=4.0 C C INPUTS ARE AS FOLLOWS C C U(1)= f (FEED)

Print file "LSTNON" Page 3 C 200 FEED=U(1) C C THE EQUATIONS ARE AS FOLLOWS.C C X(1)=WF1 C X(2)=WF2 C X(3)=WC C X(4)=WF C X(5)=F C X(6)=THETA1 C X(7)=THETA2 C XP(1)= (V/A0)*(-X(1) +Al*X(5)*COS(GAM)/(FEED*D)) XP(2)= A2*SQRT(V)*EXP(-A3/(273.+X(6))) XP(3)= A4*X(5)*V*EXP(-A5/(273.+X(7))) XP(4)= 0.0 X(4) = X(1)+X(2) XP(5)= 0.0 X(5) = (A6* (FEED**Z1) * (1.-A7*GAM)-A8-A9*V) *D+ $ A 10o*D*X (4) -Al 1 *X (3) XP(6)= 0.0 X(6) - A12*(V**Z2)*(FEED**Z3)+A13*(X(4)**Z4) XP(7)= 0.0 X(7) = A14*X(5)*(V**Z5)*(FEED**Z6)*(D**Z7) XP(8)=0.0 X(8) = A10*D*X(4)-All*X(3) C C RETURN END C C C SUBROUTINE CONTRL(T,IXC,IYC,Y, IUC,U, IREF) COMMON / INOUT / IN, I OUT COMMON/CLKBLK/HOUR, MINUTE, SECOND,MILSEC C ********** DO NOT CHANGE THESE DIMENSIONS ********* DIMENSION IXC(1),IYC(1),IUC(1),IREF(1),Y(1),U(1) C C ********** SUBSTITUTE THE VALUE OF N ************** C DIMENSION THETA(6),PHI(6),FGAIN(6,6),XBAR(3) DIMENSION TRANS(3,3),XHAT(3),XOLD(3) C C WHERE N IN THE ABOVE EQUATIONS IS 2*NEQN C INTEGER HOUR,MINUTE, SECOND,MILSEC INTEGER DCVAL C C ********** DELETE THIS SET OF CONTROL EQUNS. AND ENTER YOUR'S **** C C NOW SET UP THE MATRICES C C ********** DETERMINE THE VALUE OF N, DCVAL ********************** C N=6 DCVAL=2048+120 C 61

Print file "LSTNON" Page 4 C MINUTE=0 C SECOND=0 C MILSEC=0 C CALL STRCLK C C C SET THE THETA AND GAIN MATRICES C IF(T.GE.0.05) GOTO 230 MSTEP=0 S1=10.0 S2=1.0 DO 10 I=1,N PHI(I) =0.0 THETA (I)=(FLOAT (IXC (I) ) ) K=N+I XBAR(I)=(FLOAT (IXC (K))) 10 CONTINUE C DO 11 J=1,3 XHAT (J) =0.0 XOLD(J)=0.0 11 CONTINUE C DO 12 I=1,N DO 12 J=1,N FGAIN(I, J) =0.0 FGAIN(I, I) =1000.0 IF (I.LE.3) FGAIN(I,I)=1.0 12 CONTINUE C 230 CALL PARAM1(MSTEP,S1,S2,Y(1),U(1),)FGAIN,PHI,THETA, ERROR) C CALL STATE (MSTEP, Y (1), U ( 1 ), THETA, XBAR) C GOTO 221 IF(T.LT.0.15) GOTO 221 V=150.0 D=2.0 CALL BACK(TV, U (1),D, THETA(4), XBAR,TRANS, XHAT) C C DO 21 J=1,3 C IF (XHAT(J).L2O.O.OR.XHAT(J).GT. 0.5) XHAT (J)=XOLD (J) C XOLD (J) =XHAT(J) C 21 CONTINUE C C PRINT OUT THE VALUES C C ICOUNT=MSTEP /40 C JCOUNT=ICOUNT* 40 C WRITE (IOUT, 205) C IF (JCOUNT.EQ.MSTEP) WRITE(7,205) 205 FORMAT(//,25X,'THE TRANSFORMED STATES',//) WRITE(IOUT,206) (XHAT(J),J=1,3) C IF(JCOUNT.EQ.MSTEP) WRITE(7,206) (XHAT(J),J=1,3) 206 FORMAT (2X, 3(5X,G16.8),/) C C LET'S CALCULATE THE NEW INPUT C 221 IF(T.LT.100.0) GOTO 200 62

Print file "LSTNON" Page 5 IUC(1)=DCVAL K=IFIX(T/0.2) I=K/4 J=I*4 IF(J.EQ.K) IUC(1)=DCVAL+10 L=I/2 J=L*8 IF(J.EQ.K) IUC(1)=DCVAL-10 GOTO 350 200 IUC(1)=DCVAL 350 U(1)=(FLOAT(IUC(1) -2048))*0.0025 C C NOW UPDATE THE PHI MATRIX C YTEST=Y(1)/S1 UTEST=U(1) *S2 CALL SETUP(N,PHI,PHI,YTEST,UTEST) C C NOW BACK TRANSFER THE VALUES C DO 15 I=1,N IXC(I)=IFIX(THETA(I) *100.) 15 CONTINUE C DO 18 J=1,3 I=J+N IXC(I) =IFIX(XBAR (J)) 18 CONTINUE C C STOP THE CLOCK C C CALL STPCLK C WRITE (IOUT, 100) MINUTE, SECOND,MILSEC C100 FORMAT(' TIME IS:',I2,'.',2,'.',I5) C C PREDICT THE OUTPUT C RETURN END C SUBROUTINE YEQUNS(T,X,Y,U) C ********** DO NOT CHANGE THESE DIMENSIONS ********* DIMENSION X(1),Y(1),U(1) C ********** DELETE THIS SET OF EQUNS. AND ENTER YOUR'S **** C C C GENERATE THE NOISE C C IF(T.GT.0.05) GOTO 500 C C 11=0 C I2=0 C C 500 XNOISE=RAN(I1,I2)*50.0 C XNOISE=0.0 Y(1)=X(8)+XNOISE RETURN END 63

Print file "PARAM1" Page 1 SUBROUTINE PARAM1 (MSTEP,S1, S2,Y,U,FGAIN,PHI, THETA,ERROR) C C This is a single precision least squares parameter estimation routine C using a diagonal gain matrix and parameter constraints C COMMON/INOUT/IN, IOUT C DIMENSION THETA(6),PHI(6),FGAIN(6,6) DIMENSION BETA(6),C(6),V(6) DIMENSION PAT1 (6),PAT2 (6),PAT3 (6) C REAL* 4 LAMDA1,LAMDA2 DATA LAMDA1/ 0.95 /, LAMDA2 /1.0 / C N=6 C DO 10 I=1,3 PAT3 (I)=THETA(I) *S1 10 CONTINUE DO 11 I=4,6 PAT3(I)=THETA(I) /S2 11 CONTINUE C C LET'S CALCULATE E(K+1) C MSTEP=MSTEP+1 CALL VVSCAL(N,PAT3,PHI, Z1) ERROR=Y-Z1 C C UPDATE D,PAT3 C DO 30 J=1,N PAT1 (J)=FGAIN (J, J) *PHI (J) 30 CONTINUE C DO 31 J=1,N K=J-1 IF(K.EQ.0) GOTO 500 BETA(J) =BETA (K) +PHI (J) *PAT1 (J) GOTO 501 500 BETA(J)=LAMDA1+PHI (J) *PAT1 (J) 501 V(J)=PAT1(J) 31 CONTINUE C C Update the gain matrix C DO 32 J=1,N FGAIN (J, J) =FGAIN (J, J) -PAT1 (J) *PAT1 (J) /BETA (N) 32 CONTINUE C DO 33 I=1,N PAT2 (I)=V(I)/BETA(N) 33 CONTINUE C C UPDATE THETA C DO 13 I=1,N PAT3(I)=PAT3(I) +PAT2 (I)*ERROR 13 CONTINUE C C Impose the constraints 64

Print file "PARAMi" Rage 2 C C IF(IFLAG.EQ.0) GOTO 352 DSET1=-3.10*S1 USET1=-2.90*S1 IF(PAT3(1).LT.DSET1.OR.PAT3(1).GT.USET1) PAT3(1)=-3.0*S1 DSET2-2.90*S1 USET2=3.10*S1 IF(PAT3(2).LT.DSET2.OR.PAT3(2).GT.USET2) PAT3(2)= 3.0*S1 DSET3=-1.20*S1 USET3=-0.80*S1 IF(PAT3(3).LT.DSET3.OR.PAT3(3).GT.USET3) PAT3(3)=-1.O*S1 USET5=-1.5*PAT3 (4) DSET5=-2.5*PAT3 (4) IF(PAT3(5).LT.DSET5.OR.PAT3(5).GT.USET5) PAT3(5)=-2.0*PAT3(4) USET6=1.5*PAT3(4) DSET6=0.5*PAT3 (4) IF(PAT3(6).LT.DSET6.OR.PAT3(6).GT.USET6) PAT3(6)=PAT3(4) C C RESCALE THE THETA VECTOR C 352 DO 21 J=1,3 21 THETA(J) =PAT3 (J) /S1 DO 22 I=4,6 22 THETA(I) =PAT3 (I)*S2 C C PRINT OUT THE VALUES C ICOUNT=MSTEP / 20 JCOUNT=ICOUNT *20 WRITE(6,101) MSTEP,Y,ERROR,U IF(JCOUNT.EQ.MSTEP) WRITE(7,101) MSTEP,Y,ERROR,U 101 FORMAT(,///,' MSTEP=',I4,/,' Y=',G16.8,10X,'ERROR=',G16.8, $' U=',G16.8) WRITE(6,102) (THETA(K),K=1,6) IF(JCOUNT.EQ.MSTEP) WRITE(7,102) (THETA(K),K=1,6) 102 FORMAT(' THETA=',3(1X,G16.8),/,7X,3(1X,G16.8),//) C RETURN END 65

Print file "STATE3" Page 1 SUBROUTINE STATE (MSTEP, Y, U, THETA, XBAR) C C This a state observer with fixed gain values C COMMON/INOUT / IN, IOUT DIMENSION XBAR(3),THETA(6),VEC2(3),GAIN(3),E(3,3),F(3) C GAIN (1) =-THETA ( 1 ) GAIN (2) =-THETA(2) GAIN (3) =-THETA(3) C E (1, 1)=-THETA(1) E(1,2)=1.0 E(1,3)=0. 0 E (2,1)=-THETA(2) E(2,2)=0.0 E (2,3)=1.0 E (3,1) =-THETA(3) E(3,2)=0.0 E(3,3) =0.0 C F(1)= THETA(4) F(2)= THETA(5) F(3)= THETA(6) C ERROR=YPAST-XBAR ( 1 ) C CALL MATVEC (3,3, E, XBAR, VEC2) C DO 3 J=1,3 XBAR (J) =VEC2 (J) +F (J) *U+GAIN (J) *ERROR 3 CONTINUE C C PRINT OUT THE VALUES C ICOUNT=MSTEP / 20 JCOUNT= ICOUNT * 2 0 WRITE(IOUT, 102) ERROR, (XBAR(I), I=1, 3) IF(JCOUNT.EQ.MSTEP) WRITE (7, 102) ERROR, (XBAR(I), I=1,3) 102 FORMAT(' ERROR=',G16.8,/,' XBAR=', 3(2X,G15.8),/) YPAST=Y RETURN END 66

Print file "BACK" Page 1 SUBROUTINE BACK(TIME,V,F,D,BETA1,XBAR,T,X) C C THIS SUBPROGRAM COMPUTES THE TRANSFORMATION MATRIX C AND TRANSFORMS THE ESTIMATED STATES INTO ORIGINAL STATES C DIMENSION T(3,3),XBAR(3),X(3),B(3),C(3) C COMMON/INOUT / IN, I OUT C C NOW COMPUTE THE T MATRIX C T(1,1)= 500.0*D T(1,2)= 500.0*D T(1,3) —2000.0 C T (2,1) =-1000.0*D- (1. 783380E-06* (TIME**0. 964060) * $ (V**I1.38290)*(F**5.20270)*(D** (-1.46500))* $ (BETA1**2.60360)) T(2,2)= (-995.0+0. 050*(V-100.0)) *D-(2. 763535E-06* $ (V**2.00580)*(F**0.154940)*(D**1.00960)* $ (BETA1**(-0.0177230))) T (2,3)=(3980.0-0.20*(V-100.0))+(7.497504E-02* $ (TIME**0.322780)*(V**(0.68170))*(F**3.45220)* $ (D**(-2.02130))*(BETA1**2.12400)) C T(3, 1)=500.0*D+(2.545777E-11*(TIME**0.982830)* $ (V**3.50660)*(F**3.00210)*(D**(0.700920))* $ (BETA1**0.401560)) T (3,2) = (495.00-0.050* (V-100.0)) *D+(2. 714236E-06* $ (TIME**0.00127560)*(V**2.00880)*(F**0.154460)* $ (D**1.01130)*(BETA1**(-0.0193010))) T(3,3)=(-1980.0+0.20*(V-100.0))-(7.386619E-02* $ (TIME**0.320170)*(V**(0.681920))*(F**3.43090)* $ (D** (-2.00870)) *(BETA1**2.10970)) C C WRITE(7,201) C 201 FORMAT(//,30X,' INVERSE OF THE TRANSFORMATION MATRIX',//) C DO 10 I=1,3 C 10 WRITE(7,202) (T(I,J),J=1,3) C 202 FORMAT (3 (8X, G20.12) ) C C COMPUTE THE TRANSFORMATION MATRIX C CALL MATINV(3,3,3,T) C C WRITE(7,203) C 203 FORMAT(//,30X,'THE TRANSFORMATION MATRIX',//) C C DO 11 I=1,3 C 11 WRITE(7,204) (T(I,J),J=1,3) C 204 FORMAT(3(8X, G20.12)) C C NOW CHECK THE TRANSFORMATION MATRIX C CALL MATVEC(3,3,TXBARX) C WRITE(7,205) C 205 FORMAT(//,25X,'THE TRANSFORMED STATES',//) C WRITE(7,206) (X(J),J=1,3) C 206 FORMAT(2X,3(5X,G16.8),/) C 67

Print file "BACK" Page 2 RETURN END 68

Print file "SMAIN" Page 1 C C This program tests the on-line tool wear estimation method C on an actual system. C C C by C Kourosh Danai C INTEGER RANGE, CHANNL, GAIN, VALUE INTEGER HOUR,MINUTE, SECOND,MILSEC REAL*4 YES,ANSWER C DIMENSION THETA(6),PHI(6),FGAIN(6,6) C COMMON/ INOUT/ IN, IOUT COMMON/CLKBLK/HOUR, MINUTE, SECOND, MILSEC C DATA YES/'Y'/ DATA IN/5/,IOUT/6/ C CALL ASS IGN (3,' DLO:KDDAT.DAT' ) CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') CALL ASSIGN(7,'LP:') C CALL SINPUT (DELT, U VOLTAG, CHANNL, SCALE, S1, S2) C M=0 N=6 C C SET THE THETA AND GAIN MATRICES C DO 1 I=1,N 1 THETA(I)=0. 0 C DO 2 I=1,N DO 2 J=1,N FGAIN(I, J) =0.0 FGAIN (I, I) =1000.0 IF(I.LE.3) FGAIN(I, I)=1.0 2 CONTINUE C C determine the gain for the ADC board C CALL KDADC (VOLTAG, CHANNL, GAIN) C GOTO 10 20 CALL STPCLK WRITE(IOUT,101) M 101 FORMAT(' Would you like to CONTINUE the program? Y/N?!', $' MSTEP=',13,/,) READ(IN, 102) ANSW 102 FORMAT (A2) IF(ANSW.NE.YES) GOTO 40 C 10 CALL DRVRED ( 1,VALUE) IF(VALUE.EQ.-2) GOTO 10 C C Set the flag for the constraints C 69

Print file "SMAIN" Page 2 IFLAG=1 C C Reset the PHI vector and FGAIN matrix C DO 3 I=1,N PHI(I)=0.0 3 CONTINUE C YTEST=0.0 UTEST=U*S2 CALL SETUP (N, PHI,PHI, YTEST,UTEST) C DO 4 I=1,N DO 4 J=1,N FGAIN(I, J) =0.0 FGAIN (I, I) =1000.0 IF(I.LE.3) FGAIN(I, I)=1.0 4 CONTINUE C CALL SETTIM(0,0,0,0) CALL ADC (GAIN, IVALUE) CALL FACTOR ( SCALE, VOLTAG, IVALUE, OUTPUT) WRITE (IOUT, 120) OUTPUT 120 FORMAT(' FORCE=',G15.8) C DCVAL=OUTPUT CALL SETTIM(0,0,0,0) C CALL STRCLK C 30 IF(SECOND.LT.DELT) GOTO 30 CALL SETTIM(0,0,0,0) CALL ADC (GAIN, IVALUE) CALL FACTOR(SCALE, VOLTAG, IVALUE, OUTPUT) Y-OUTPUT-DCVAL C CALL PARAM1i (M, S1, S2,Y,U, FGAIN, PHI, THETA, ERROR, IFLAG) D IFLAG=0 C C Store the variables C WRITE(3) M,Y,U, (THETA(J),J=1, 6),ERROR C C Now update the PHI vector C YTEST=Y/S1 UTEST=U*S2 CALL SETUP (N, PHI,PHI, YTEST,UTEST) C WRITE (IOUT, 100) MINUTE, SECOND,MILSEC 100 FORMAT(' TIME IS:' I2'' I2'.',I5) C CALL DRVRED ( 1, VALUE).IF(VALUE.EQ.-2) GOTO 20 GOTO 30 C STOP THE CLOCK 40 WRITE(6,103) 103 FORMAT(' Would you like to print out the values?!, Y/N?',/) READ(5,102) ANSW 70

Print file "SMAIN" Page 3 IF(ANSW.NE.YES) GOTO 70 C WRITE (7,104) 104 FORMAT(1X,'Mstep',7X,' Output', 9X,'Input', /,'Alphal', 9X, $'Alpha2',9X,'Alpha3',9X,'Betal', 10X,'Beta2',10X,'Beta3',10X, $'ERROR' //!) C REWIND 3 DO 5 J=1,M READ(3) MSTEP,Y,U, (THETA(I),I=1, 6),ERROR WRITE(7,105) MSTEP,Y,U, (THETA(I),I=1, 6),ERROR 105 FORMAT(1X, I4,7X,2(1X, G14.8),/,7(1X, G14.8)) 5 CONTINUE C 70 CLOSE (UNIT=3) STOP END 71

Print file "SAMPLE" Page 1 C C This program samples the forces one axis at a time C C C by C Kourosh Danai C INTEGER RANGE, CHANNL, VALUE, DELT INTEGER GAIN1, GAIN2, GAIN3 INTEGER CH1, CH2, CH3 INTEGER HOUR,MINUTE,SECOND,MILSEC REAL*4 YES,ANSWER C COMMON/INOUT/IN, IOUT COMMON/CLKBLK/HOUR,MINUTE,SECOND,MILSEC C DATA YES/'Y' / DATA IN/5/,IOUT/6/ C CALL ASSIGN(1,'DL0:KDX1.DAT' ) D CALL ASSIGN(2,'DL0:Y1.DAT' ) D CALL ASSIGN(3,'DL0:Z1.DAT' ) CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') CALL ASSIGN(7,'LP:') C CALL SINPUT(DELT,U,VOLTAG,CHANNL,SCALE,S1, S2) C WRITE(IOUT,505) DELT 505 FORMAT(' DELT=',I3) M=0 CH1=1 D CH2=2 D CH3=3 C C determine the gain for the ADC board C CALL KDADC(VOLTAG, CHl1, GAIN1) D CALL KDADC(VOLTAG,CH2,GAIN2) D CALL KDADC(VOLTAG,CH3, GAIN3) C GOTO 10 20 CALL STPCLK WRITE(IOUT, 101) M 101 FORMAT(' Would you like to CONTINUE the program? Y/N?!', $' MSTEP=',I3,/,) READ (IN, 102) ANSW 102 FORMAT(A2) IF(ANSW.NE.YES) GOTO 40 C 10 CALL DRVRED ( 1, VALUE) IF(VALUE.EQ.-2) GOTO 10 C CALL SETTIM(0,0,0,0) CALL ADC(GAIN1,IVAL1) D CALL ADC(GAIN2,IVAL2) D CALL ADC(GAIN3,IVAL3) CALL FACTOR(SCALE, VOLTAG, IVAL1,OUT 1) D CALL FACTOR(SCALE,VOLTAG, IVAL2,OUT2) 72

Print file "SAMPLE" Page 2 D CALL FACTOR (SCALE, VOLTAG, IVAL3, OUT3) C WRITE(IOUT, 120) OUT1 D WRITE (IOUT,120) OUT1,OUT2,OUT3 120 FORMAT(' XFORCE=',G16. 8//) D 120 FORMAT(' XFORCE=',G15.8,2X,'YFORCE=',G16.8,2X,'ZFORCE=',G16.8//) C DCVAL1=OUT1 D DCVAL2=OUT2 D DCVAL3=OUT3 C CALL SETTIM(0,0,0,0) C CALL STRCLK C 30 IF(MILSEC.LT.DELT) GOTO 30 CALL SETTIM(0,0,0,0) C CALL ADC(GAIN2,IVAL2) D CALL ADC(GAIN2,IVAL2) D CALL ADC(GAIN3,IVAL3) C CALL FACTOR(SCALE,VOLTAG, IVAL1,OUT1) D CALL FACTOR(SCALE, VOLTAG, IVAL2,0 OUT2) D CALL FACTOR (SCALE, VOLTAG, IVAL3, OUT3) C D X=OUT1 X=OUT1-DCVAL1 D Y=OUT2-DCVAL2 D Z=OUT3-DCVAL3 C C Store the variables C M=M+1 WRITE(1) M,X D WRITE(2) M,Y D WRITE(3) M,Z C WRITE(IOUT,107) M, X 107 FORMAT(' M=',I3,5X,' FORCE=',G16.8) C D WRITE (IOUT, 100) MINUTE, SECOND,MILSEC D 100 FORMAT(' TIME IS:',I2,'.',I2,'',I5) C CALL DRVRED(1,VALUE) IF(VALUE.EQ.-2) GOTO 20 GOTO 30 C C STOP THE CLOCK C 40 WRITE(6,103) 103 FORMAT(' Would you like to print out the values?!, Y/N?',/) READ(5,102) ANSW IF (ANSW.NE.YES) GOTO 70 WRITE(7,104) 104 FORMAT(1X,'Mstep',7X,'XFORCE'//) D104 FORMAT(1X,'Mstep',7X,'XFORCE' 9X,'YFORCE' 9X,' ZFORCE'//) C 73

Print file "SAMPLE" Page 3 REWIND 1 D REWIND 2 D REWIND 3 C DO 5 J=1,M READ(1) MSTEP,X D READ(2) MSTEP,Y D READ(3) MSTEP,Z WRITE(7,105) MSTEP,X D WRITE(7,105) MSTEP,X,Y, Z 105 FORMAT(1X, I4,7X, Gl4. 8) D 105 FORMAT(1X, I4,7X, 3(lX, G14.8)) 5 CONTINUE C 70 CLOSE (UNIT=l) D CLOSE (UNIT=2) D CLOSE (UNIT=3) C STOP END 74

Print file "ESTIM" Page 1 C C This program tests the on-line tool wear estimation method C on the collected data. C C C by C Kourosh Danai C REAL* 4 YES,ANSWER C DIMENSION THETA(2),PHI(2), FGAIN(2,2) C COMMON/INOUT/IN, IOUT C DATA YES/'Y' / DATA IN/5/,IOUT/6/ C CALL AS S IGN ( 1,' OPARX2.DAT' ) CALL ASSIGN (2,' KDX2.DAT' ) C CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') CALL ASSIGN(7,' LP:' ) C S1=10.0 S2=1.0 U =.254 MVAL= 8 4 N=2 E=0.0 C C SET THE THETA VECTOR AND GAIN MATRIX C DO 1 I=1,N THETA(I)=. 0 PHI(I)=0.0 1 CONTINUE C DO 2 I=1,N DO 2 J=1,N FGAIN(I,J)=0.0 FGAIN (I, I) =100.0.0 IF(I.EQ.1) FGAIN(I,I)=1.0 2 CONTINUE C REWIND 2 C DO 11 K=1,MVAL READ(2) MSTEP, Y IF(MSTEP.GE.17) Y = Y + 21.88 IF(MSTEP.GE.33) Y = Y + 55.88 IF(MSTEP.GE.50) Y = Y + 21.49 IF(MSTEP.GE.63) Y = Y + 13.28 IF(MSTEP.GE.50) Y = Y + 23.05 C WRITE (IOUT, 120) Y C 120 FORMAT(' FORCE=',G15.8) C CALL OPAR21 (MSTEP, S1, S2,Y, U, FGAIN, PHI, THETA, ERROR, E) D IFLAG=0 CALL STATE1 (Y, U, THETA, XBAR) X = XBAR/200.0 75

Print file "ESTIM" Page 2 WRITE(IOUT,109) X 109 FORMAT(' The estimated state=',G16.8,/) C C Store the variables C WRITE(1) MSTEP, Y, U, (THETA(J), J=1,2),ERROR, X C C Now update the PHI vector C YTEST=Y/S 1 UTEST=U*S2 CALL SETUP (N, PHI,PHI, YTEST,UTEST) 11 CONTINUE C CLOSE (UNIT=1) CLOSE (UNIT=2) C STOP END 76

Print file "SINPUT" Page 1 SUBROUTINE SINPUT (DELT,U,VOLTAG,CHANNL, SCALE, Sl, S2) C C Subroutine SINPUT is designed to display various sampling parameters for the C user, and to prompt him for new ones. The parameters are passed back to the C calling program through the common blocks. C INTEGER DELT, CHANNL, RANGE C COMMON / INOUT / IN, IOUT C C Increment the parameter flag. The parameter flag is used to tell the proC gram whether this is the user's first run through the program. If it is C the first time through the program, all the sampling parameters will be inC itialized by the subroutine SINPUT. On subsequent passes through the proC gram, the sampling parameters that the user entered during the previous pass C will be retained. Thus it will be easier for the user if he wants to change C only one parameter per pass. C PARFLG = PARFLG + 1 C C Check to see if this is the first time the user is passing through SINPUT, C and branch appropriately. C IF ( PARFLG.GT. 1 ) GO TO 9 C C Since this is the user's first pass through the program, initialize cutC ting parameters. Note that the ADCGAN factor has units of volts/unit number. C DELT = 500 U = 0.01 RANGE = 4 CHANNL = 1 SCALE = 1250.0 S1=10.0 S2=1.0 U=U*25.4 C C C Display a header. C 9 IF(RANGE.EQ.1) VOLTAG=10.24 IF(RANGE.EQ.2) VOLTAG=5.12 IF(RANGE.EQ.3) VOLTAG=2.56 IF(RANGE.EQ.4) VOLTAG=1.28 WRITE(IOUT, 101) 101 FORMAT('1'/'1'/'1'/ + 20X,'TOOL WEAR ESTIMATION PROGRAM' C C Display the present value of all the sampling parameters. C WRITE(IOUT, 102) DELT, U, VOLTAG, CHANNL, SCALE, S1, S2 102 FORMAT( ////// 16X,'SAMPLING PARAMETERS DISPLAY:' // + 2X,' 1) Sampling period is:', + I7,' (milsec)' / + 2X,' 2) Feed input is:', + G10.4,' (rmm)' / + 2X,' 3) Output range of the ADC channel is', + G10. 4,' (volts)' / + 2X,' 4) Channel no. on the ADC board is:', + I7,' (unitless)' / + 2X,' 5) Scale factor of the charge amplifier is:' 77

Print file "SINPUT" Page 2 + G15.8,' (Newtons/Volt)' / + 2X,' 6) Force dividing factor for estimation: + G7.2,' (unitless)' / + 2X,' 7) Input multiplication factor for estimation is', + G7.2,' (seconds)' /////) C C Ask the user if he wants to change any of the present values of the C sampling parameters. C WRITE(IOUT, 103) 103 FORMAT(' To change any of the parameter values, enter the', +' number of the parameter.'/ +' Otherwise RETURN to continue READ(IN, 104, ERR=9) IFLAG 104 FORMAT (I2) C IF ( IFLAG.LT. 0.OR. IFLAG.GT. 7 ) GO TO 9 C IF ( IFLAG.EQ. 0 ) GO TO 5000 C GO TO ( 301, 302, 303, 304, 305, 306, 307 ) IFLAG C 301 ASSIGN 301 TO I WRITE ( IOUT, 201) 201 FORMAT ('1'/'1'/'1'/'$Enter a new value for the', +' sampling period:' ) READ ( IN, 202, ERR=990 ) DELT 202 FORMAT ( I2 ) GO TO 9 C 302 ASSIGN 302 TO I WRITE ( IOUT, 203) 203 FORMAT ('1'/'1'/'1'/'$Enter a new value for the', +' feed input (in/rev):' ) READ ( IN, 204, ERR=990 ) U 204 FORMAT(G8.4) U=U*25.4 GO TO 9 C 303 ASSIGN 303 TO I WRITE ( IOUT, 205) 205 FORMAT ('1'/'1'/'1'/' Enter a new value for the', +' voltage range on the ADC board:' / +' 1. +/- 10.24 V',/, +' 2. 5.12 V',/, +' 3. +/- 2.56 V',/, +' 4. +/- 1.28 V',/, +' Please enter 1, 2, 3, or 4',/) READ ( IN, 202, ERR=990 ) RANGE GO TO 9 C 304 ASSIGN 304 TO I WRITE ( IOUT, 206 ) 206 FORMAT ('1'/'1'/'1'/'$Enter a new value for the', +' channel No.:' ) READ ( IN, 202, ERR=990 ) CHANNL GO TO 9 C 305 ASSIGN 305 TO I WRITE ( IOUT, 207 ) 207 FORMAT ('1'/'1'/'1'/'$Enter a new value for the', 78

Print file "SINPUT" Page 3 +' amplifier scale factor:' READ ( IN, 204, ERR=990 ) SCALE GO TO 9 C 306 ASSIGN 306 TO I WRITE ( IOUT, 208 ) 208 FORMAT ('1' /'1' /'1' /'$Enter a new value for the', +' dividing factor:' ) READ ( IN; 204, ERR=990 ) S1 GO TO 9 C 307 ASSIGN 307 TO I WRITE ( IOUT, 209) 209 FORMAT ('1'/'1'/'1'/'$Enter a new value for the', +' multiplication factor:' ) READ ( IN, 204, ERR=990 ) S2 GO TO 9 C C 990 WRITE(IOUT, 991) 991 FORMAT(' ***** ERROR ***** Invalid entry, try again', /,'') GO TO I C C Return to the calling program. C 5000 RETURN END 79

Print file "OPAR21" Page 1 SUBROUTINE OPAR21 (MSTEP, S 1, S2, Y, U, FGAIN, PHI, THETA, ERROR, E ) C C General Description: This is an output error estimation C routine for a first order model C C by C Kourosh Danai C C COMMON / INOUT / IN, IOUT C DIMENSION THETA(2),PHI(2),FGAIN(2,2) DIMENSION PAT1(2),PAT2(2),PAT3(2),PAT4(2),PAT5(2,2) C REAL*4 LAMDA1, LAMDA2 DATA LAMDA 1. /, LAMDA2 /1. / DATA C1/1.0/ C N=2 C PAT3(1) = THETA(1)*S1 PAT3(2) = THETA(2)/S2 C C E(1)=0.O C E(2)=0.0 C CALL VVSCAL (N, PAT3,PHI, Z1) CALL VECMAT (NNPHI, FGAIN, PAT1) CALL VVSCAL(N, PAT1,PHI, Z2) EDOT=Y-Z 1 C ERROR= EDOT+C 1 *E ERROR= ERROR/(1.+Z2) C C UPDATE THETA C DO 13 I=1,N CALL MATVEC (N, N, FGAIN, PHI, PAT2) 13 PAT3 (I)=PAT3 (I) +PAT2 (I) *ERROR C C COMPUTE THE PRIOR ERROR TERMS C CALL VVSCAL(N, PAT3,PHI, Z3) E=Y-Z3 C C UPDATE FGAIN C CALL VVMAT (N,N,PAT2,PAT1,PAT5) DO 14 I=1,N DO 14 J=1,N FGAIN (I, J) = (FGAIN (I,J)-PAT5 (I,J) / (LAMDA1/LAMDA2+Z2 ) ) /LAMDA1 14 CONTINUE C C NOW BACK TRANSFER THE VALUES THETA(1) =PAT3 (1) /S1 THETA(2) =PAT3 (2) *S2 C C Print the stuff out C WRITE(IOUT, 102) MSTEP,Y, ERROR, PHI (1),PHI (2) 80

Print file "OPAR21" Page 2 102 FORMAT(' M=',I4,/,' Y=',G16.8,' ERROR=',G16.8,/, $' PHI =',2(2X,G16.8),/) WRITE(IOUT, 103) THETA(1),THETA(2) 103 FORMAT(2X,' THETA =',2(2X, G16.8),//) C RETURN END 81

Print file "STATEl" Page 1 SUBROUTINE STATE1(Y,U,THETA, XBAR) C C This a state observer with fixed gain values C COMMON/INOUT/IN, IOUT DIMENSION THETA (2) C GAIN=-THETA ( 1 ) C E=-THETA (1) C F= THETA(2) C ERROR=YPAST-XBAR C VEC2 = E*XBAR C XBAR = VEC2 + F * U + GAIN*ERROR C WRITE(6,102) ERROR, XBAR D WRITE(7,102) ERROR, XBAR 102 FORMAT(' ERROR=',G16.8,/,' XBAR=',2X,G15.8,/) YPAST=Y RETURN END 82

Print file "KDADC" Page 1 C ANALOG TO DIGITAL SOURCE CODE INTEGER IADC, ICH, IGAIN C REAL*4 Y S1 DATA IN/5/,IOUT/6/ DATA Y/'Y' / C CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') C 155 WRITE(IOUT, 110) 110 FORMAT(3X,'WOULD YOU LIKE TO READ THE I/O PORT', /, $' REGISTER, Y/N?!',/) READ(IN, 123) S1 123 FORMAT (A2) C IF(S1.NE.Y) GOTO 10 C CALL DRVRED ( 1, NUMBER) WRITE (IOUT, 124) NUMBER 124 FORMAT(3X,' THE I/O PORT NO. IS =', I6, /) GOTO 155 C 10 WRITE(IOUT, 111) 111 FORMAT(3X,'>>>ANALOG RANGE:', /, Z3X,'1. +/-10 V' /, Z3X,'2. +/-5 V', / Z3X,'3. +/-2.5 V' /, Z3X,'4. +/-1.25 V',/, Z3X,'>>>PLEASE ENTER 1,2,3, OR 4?') READ (IN, 211) IC 211 FORMAT(I2) GOTO(310,320,330,340),IC GOTO 10 310 IGAIN=0 GOTO 20 320 IGAIN=4 GOTO 20 330 IGAIN=8 GOTO 20 340 IGAIN=12 C 20 WRITE (IOUT, 121) 121 FORMAT(3X,'~>>>PLEASE SPECIFY THE CHANNEL #?0-7' ) READ(IN,211)IC IF(IC.EQ. 0)IC=8 GOTO(410,420,430,440,450,460,470,480),IC GOTO 20 410 ICH=2**8 GOTO 30 420 ICH=2**9 GOTO 30 430 ICH=2**9+2**8 GOTO 30 440 ICH=2**10 GOTO 30 450 ICH=2**10+2**8 GOTO 30 460 ICH=2**10+2**9 GOTO 30 470 ICH=2**10+2**9+2**8 8 3

Print file "KDADC" Page 2 GOTO 30 480 ICH=2**11 C 30 IVAL=ICH+IGAIN 40 CALL ADC (IVAL, IADC) WRITE (IOUT, 131) IADC 131 FORMAT (3X,'ANALOG TO DIGITAL VALUE=',I8,/, Z3X,'>>WOULD YOU LIKE TO:',/, Z3X,' 1. CALL ADC AGAIN;',/, Z3X,' 2. REDEFINE CHANNEL AND ANALOG RANGE;',/, Z3X,' 3. QUIT...?',/, Z3X,'...PLEASE ENTER 1,2,3...') READ(IN,211)IC GOTO(40,10,50),IC 50 CONTINUE STOP END 84

Print file "ADC" Page 1.TITLE ADC.GLOBL ADC ADC: MOV (R5)+,RO MOV @ (R5) +, @#170400 INC @#170400 3$: TSTB @#170400 BPL 3$ MOV @#170402,@ (R5)+ RTS PC.END 85

Print file "SSUBS" Page 1 SUBROUTINE KDADC (VOLTAG, CHANNL, GAIN) C C This program computes the gain for the ADC.MAC C INTEGER RANGE,CHANNL,GAIN C COMMON /INOUT/IN,IOUT C IF(VOLTAG.EQ. 1.28) IRANGE=12 IF(VOLTAG.EQ.2.56) IRANGE=8 IF(VOLTAG.EQ. 5.12) IRANGE=4 IF(VOLTAG.EQ.10.24) IRANGE=0 C IF(CHANNL.EQ.0) IC=2**11 IF(CHANNL.EQ.1) IC=2**8 IF(CHANNL.EQ.2) IC=2**9 IF(CHANNL.EQ.3) IC=2**9+2**8 IF(CHANNL.EQ.4) IC=2**10 IF(CHANNL.EQ.5) IC=2**10+2**8 IF(CHANNL.EQ.6) IC=2**10+2**9 IF(CHANNL.EQ.7) IC=2**10+2**9+2**8 C GAIN=IC+IRANGE C RETURN END C C C C SUBROUTINE FACTOR (SCALE,VOLTAG,VALUE,OUTPUT) C C This subroutine computes the real force C INTEGER VALUE C D VOLT= (FLOAT (VALUE-2048)) *VOLTAG/2048 VOLT= (FLOAT (VALUE)) *VOLTAG/4095 OUTPUT=VOLT*SCALE C RETURN END C C C C SUBROUTINE SETTIM (HVAL,MINVAL, SECVAL, MILVAL) C C This subroutine sets the clock C INTEGER HVAL,MINVAL, SECVAL,MILVAL INTEGER HOUR,MINUTE, SECOND,MILSEC C COMMON/CLKBLK/HOUR, MINUTE, SECOND, MILSEC HOUR-=HVAL MINUTE=MINVAL SECOND=SECVAL MILSEC=MILVAL RETURN

Print file "SSUBS" Page 2 END SUBROUTINE MATMUL (MA, NA, MB, NB, A, B, AB) C C C THIS SUBPROGRAM MULTIPLIES THE TWO MATRICES AND PUTS THE C RESULTS IN MATRIX AB C C DIMENSION A(MA, NA),B (MB, NB),AB (MA, NB) C DO 1 I=1,MA DO 1 J=1,NB AB(I, J) =O.O DO 1 K=1,NA 1 AB(I, J) =AB(I, J) +A(I, K) *B (K,J) RETURN END C C C C SUBROUTINE MATVEC (MA, NA, A, V, W) C C C THIS SUBPROGRAM MULTIPLIES THE MATRIX A WITH THE VECTOR C V AND RETURNS THE RESULT IN VECTOR W C C C DIMENSION A(MA,NA),V(NA),W(MA) C DO 2 I=1,MA 2 W(I)=O.O DO 3 I=1,MA DO 3 J=1,NA W(I) =(I) +A(I, J) *V(J) 3 CONTINUE RETURN END C C C C C C SUBROUTINE VECMAT (MA, NA, V, A, W) C C C THIS SUBPROGRAM MULTIPLIES THE VECTOR V BY MATRIX A AND C PUTS THE RESULT IN W DIMENSION V (MA), A (MA, NA), W (NA) DO 3 J=1,NA 87

Print file "SSUBS" Page 3 W(J)=0.0 DO 3 I=1,MA W(J)=W(J)+V(I) *A(I,J) 3 CONTINUE RETURN END C C C C SUBROUTINE VVSCAL (NV,W, Z) C C C THIS SUBROUTINE MULTIPLIES THE TWO VECTORS AND PUTS THE C RESULTS IN THE SCALAR Z C C DIMENSION V(N), W (N) C Z=0.0 DO 4 I=1,N Z=Z+v(I) *W(I) 4 CONTINUE RETURN END C C C C SUBROUTINE VVMAT (NV, NW, V, W, A) C C C THIS SUBROUTINE MULTIPLIES THE TWO VECTORS V AND W C AND PUTS THE RESULT IN MATRIX A C C DIMENSION V (NV),W (NW) A (NV, NW) C C DO 6 I=1,NV DO 6 J=1,NW A(I, J)-=V(I) *W (J) 6 CONTINUE RETURN END C C C C SUBROUTINE SETUP (NOLD,FRESH, YU) DIMENSION OLD (N), FRESH (N) C I=N L-N/2+1 50 J=I-1 FRESH(I)=OLD(J) IF(I.EQ.L) FRESH(I)=U I=I-1 IF(I.EQ. 1) GOTO 60 GOTO 50 88

Print file "SSUBS" Page 4 60 FRESH(I) =-Y RETURN END C C C C C SUBROUTINE MATINV (MA, NA, N, A) C C C C Replaces matrix A by its inverse C DIMENSION A(MA,NA) C DO 1 I=1,N DUM=A(I, I) DO 3 J=-1,N 3 A(I, J)=A(I, J)/DUM A(I, I) =1. 0/DUM DO 4 J=-1,N IF(I.EQ.J) GOTO 4 DUM=A(J, I) DO 5 K=1,N 5 A(J,K)=A(J,K)-A(I,K) *DUM A(J, I)=-A(I, I) *DUM 4 CONTINUE 1 CONTINUE C RETURN END C C C C C C FUNCTION DGAMMA(N) C C C THIS FUNCTION CALCULATES THE FACTORIAL OF N C C DGAMMA=1.DO DO 1 I=1,N DGAMMA=DGAMMA* I 1 CONTINUE C RETURN END C C C SUBROUTINE TRANS(MA,NA,A,AT) C This subroutine transforms A and puts the result C in AT C C 89

Print file "SSUBS" Page 5 DIMENSION A(MA,NA), AT (NA, MA) MAT=NA NAT=MA C DO 2 I=1,MAT DO 2 J=1,NAT AT(I, J) =O. O 2 AT (I, J) =AT (I, J) +A(J, I) C RETURN END 90

Appendix C This appendix presents a set of routines to plot and print data stored in a data file. These routines are based on a plotting package developed by Leal Lauderbaugh [1] for a Digital Equipment Corporation (DEC) LSI-11/23 plus microcomputer. The components of this package are: DISPLY: Main program to call different subroutines to plot and print the stored data. PLOT: Subroutine to plot the selected variable on a Tecktronix 4006-1 graghics terminal. HRDCPY: Subroutine to generate a hard copy of the plot on an xy recorder. PRNT: Subroutine to-print out the data on either a printer or video terminal (should be modified for different data files). CHOICE: Subrutine to select the data for plotting (should be modified for different data files).

Print file "DISPLY" Page 1 C General Description: C C This package is used to print and plot different variables C C External References: C C C PLOT: routine to plot the sampled forces on the graphics terminal C CHOICE: routine to select the variable to be plotted C HRDCPY: routine to plot the sampled forces on the hardcopy unit C PRNT: routine to print out the sampled forces C INTEGER OPTION, POINTS C DIMENSION PARAM(1000) COMMON /INOUT/ IN, IOUT COMMON /ACQBLK/ POINTS,PARAM DATA IN /5/, IOUT /6/, PARAM/1000*0.0/ C C Assign devices C CALL ASSIGN(5,'TT:') CALL ASSIGN(6,'TT:') C C WRITE (IOUT, 321) 321 FORMAT(' Enter the no. of points',/) READ (IN, 322) POINTS 322 FORMAT (I3) C C Display the main menu. C 10 WRITE ( IOUT, 20 ) 20 FORMAT ('1', 35X,'MAIN MENU' ////// + 22X,'1) Plot the values on the terminal' / + 22X,'2) Print out the forces' / + 22X,'3) Quit' ///////) C C Prompt user for menu option. C WRITE ( IOUT, 30 ) 30 FORMAT ('$Choose option ( 1, 2, 3 )' ) READ ( IN, 40, ERR = 10 ) OPTION 40 FORMAT ( I2 ) C C Check to see that the user chose a valid option. C IF ( OPTION.GT. 3.OR. OPTION.LT. 1 ) GO TO 10 C C Branch to carry out the option requested C GO TO (1000,2000,3000) OPTION C C Plot the desired force. C 1000 CALL PLOT GO TO 10 C C Print out the forces. 92

Print file "DISPLY" Page 2 C 2000 CALL PRNT GO TO 10 C C Leave the program. C 3000 STOP END 93

Print file "PLOT" Page 1 SUBROUTINE PLOT C C General Description: C C This subroutine plots out the desired parameter. C The user chooses which parameter he wants to print out via the plotC ting menu. PLTFRC uses the LEALPT plotting package. The plot is sent to C the Tektronix 4006-1 graphics terminal. C COMMON /INOUT/ IN, IOUT C COMMON /MXMN/ DVMAX, DVMIN, RIDVMX, RIDVMN C C DVMAX: REAL. The maximum value of the dependent variable. C DVMIN: REAL. The minimum value of the dependent variable. C RIDVMX: REAL. The maximum value of the independent variable. C RIDVMN: REAL. The minimum value of the independent variable. C COMMON /RESBLK/ DV, RIDV, LOWPT, HIGHPT C C DV: REAL. An array of dimension 1000. Stores the deC pendent variable array for use in the plotC ting package. C RIDV: REAL. An array of dimension 1000. Stores the inC dependent variable array for use in the plotC ting package. In this subroutine, a plotC ting index, from 0 up to the number of C points plotted, in increments of one, is C stored as the dependent variable. C LOWPT: INTEGER. The lowest data point to plot, in terms of C the numerical value of the independent C variable RIDV. LOWPT and HIGHPT are used C for "windowing" operations with the plotting C routines, i.e. the user can look at all of C his data first, then pick a smaller portion C of the independent axis to look at. C HIGHPT: INTEGER. The highest data point to plot, in terms of C the numerical value of the independent C variable RIDV. LOWPT and HIGHPT are used C for "windowing" operations with the plotting C routines, i.e. the user can look at all of C his data first, then pick a smaller portion C of the independent axis to look at. C C C External References: C C C AXES: routine to draw the axes of the plot C AXTYPE: routine to choose which axes to draw, i.e. which quadrants C DRAW: routine to draw a line on the plot from the current position C to the coordinates specified as arguments C MAXMIN: routine to find the maximum and minimum values of the C dependent and independent variable arrays C MOVE: routine to move from the current screen position to the posiC tion specified in the routine arguments C PAGE: routine to clear the screen of the graphics terminal C SCALE: routine to calculate the scale factors and offsets for the C plotting routines C C 94

Print file "PLOT" Page 2 C C C Declaration statements. C COMMON /ACQBLK/ POINTS,PARAM C INTEGER LOWPT, HIGHPT, OPTION, POINTS LOGICAL*1 ANSWER, Y C DATA Y/'Y'/ C C Dimension statements. C DIMENSION DV ( 1000 ), RIDV ( 1000 ), PARAM(1000) C C Initialize the plotting parameters. C LOWPT = 1 HIGHPT = POINTS C C Display a header. C 10 WRITE (IOUT, 12) 12 FORMAT (' 1'/'1'/'1' / + 29X,'PLOTTING PROGRAM' C C Display the plotting menu. C WRITE ( IOUT, 20 ) 20 FORMAT('1', 32X,'PLOTTING MENU' /// + 12X,'1) Select a new variable to Plot on the graphics terminal', +,/, + 12X,'2) Repeat the previus plot',/, + 12X,'3) Clear the screen of the graphics terminal',/, + 12X,'4) Change the plotting window',/, + 12X,'5) Return to the main menu' ///// ) C C Prompt user for menu option. C WRITE ( IOUT, 30) 30 FORMAT ('$Please choose option ( 1, 2, 3, 4, 5 ) READ ( IN, 40, ERR = 10 ) OPTION 40 FORMAT ( I2 ) C C Check to see that the user chose a valid option. C IF ( OPTION.GT. 5.OR. OPTION.LT. 1 ) GO TO 10 C C Branch to carry out the option requested C GO TO ( 1000, 2000, 3000, 4000, 5000 ) OPTION C C Put the force numbers that the user wants to print out into a one C dimensional array. C 1000 CALL CHOICE LOWPT=1 HIGHPT = POINTS C 444 DO 100 I = LOWPT, HIGHPT C 95

Print file "PLOT" Page 3 DV ( (I - LOWPT) + 1 ) = PARAM(I) RIDV ( (I- LOWPT) + 1 ) = (I - LOWPT) C 100 CONTINUE C C Find the maximum and minimum values of the plotting arrays DV, and RIDV C for use in scaling the graphs. C 1010 CALL MAXMIN ( (HIGHPT - LOWPT), DV, RIDV ) C C Find the axes type. C 1020 CALL AXTYPE ( IQUAD) C C Find the scale for the plot. C CALL SCALE( IQUAD, DVSCL, IDVOFF, RIDVSC, IIDVOF ) C C C See if the scale need be changed C WRITE ( IOUT, 200 ) DVMAX, DVMIN, RIDVMX, RIDVMN 200 FORMAT ('1'/'1'/'1'/ +' The maximum value of the dependent variable is', G11.4, /, +' The minimum value of the dependent variable is', G11.4, /, +' The maximum value of the indepent variable is',G11.4,/, +' The minimum value of the indepent variable is',G11.4,/, +,/,' Do you wish to make any changes? Y/N?') C READ(IN, 101) ANSWER 101 FORMAT (A2) IF(ANSWER.EQ.Y) GOTO 329 GOTO 333 C 329 WRITE (IOUT, 520) 520 FORMAT(' Input the new DVMAX,DVMIN, RIDVMX,RIDVMN' ) READ (IN, 530) DVMAX, DVMIN, RIDVMX, RIDVMN 530 FORMAT(4(G10.4)) GOTO 1020 C C Draw the axes. C 333 CALL AXES ( IQUAD) C C Plot the forces. C 555 DO 400 K = LOWPT, HIGHPT C IX = IFIX ( RIDVSC * RIDV ( (K - LOWP.T) + 1 ) ) + IIDVOF IY = IFIX ( DVSCL * DV ( (K - LOWPT) + 1)) + IDVOFF C IF ( K.NE. LOWPT ) GO TO 390 CALL MOVE ( IX, IY ) GO TO 400 C 390 CALL DRAW ( IX, IY ) C 400 CONTINUE C 4099 WRITE ( IOUT, 4100 ) 4100 FORMAT ( /'1'/'1'/'1'/'$Press RETURN to continue:' ) 96

file "PLOT" Page 4 READ ( IN, 4200 ) ANSWER 4200 FORMAT ( Al ) C C See if a hardcopy of the plot is needed. C WRITE (IOUT,203) 203 FORMAT( /'1'/'1'/'1'/' Would you like to get a hardcopy of' $' this plot?! Y/N?!', /) READ (IN, 101) ANSWER IF(ANSWER.EQ.Y) CALL HRDCPY GO TO 10 C C Repeat the plot C 2000 GOTO 333 C C Clear the plotting screen. C 3000 CALL PAGE GO TO 10 C C Prompt the user for new values of LOWPT and HIGHPT, C which define the plotting window. C 4000 WRITE ( IOUT, 5100 ) POINTS, LOWPT, HIGHPT 5100 FORMAT ('1'/'1'/'1'/ +' The original data sample ranges from point 1 to point', + I5,'.' // +' The plotting window currently extends from data', +' point', I5,' to data point', I5,'.' /'1'/'1'/// +'$Please enter the lowest data point to plot:' ) C READ ( IN, 5200 ) LOWPT 5200 FORMAT ( I5) C WRITE ( IOUT, 5300 ) 5300 FORMAT ( //'$Please enter the highest data point to plot:' ) C READ ( IN, 5200 ) HIGHPT C C Replot the variable C GO TO 444 C C Return to the calling program. C 5000 RETURN END 9 7

Print file "HRDCPY" Page 1 SUBROUTINE HRDCPY C C General Description: C C This subroutine gives a hard copy of the plot generated by subroutine PLOT C HRDCPY uses the LEALPT plotting package. The plot is sent to the C XY Recorder ( Hardcopy unit ). C COMMON /INOUT/ IN, IOUT COMMON /MXMN/ DVMAX, DVMIN, RIDVMX, RIDVMN COMMON /RESBLK/ DV, RIDV, LOWPT, HIGHPT C INTEGER LOWPT, HIGHPT C DIMENSION DV(1000), RIDV(1000), PARAM(1000) C COMMON /ACQBLK/ POINTS,PARAM C C Initialize the hardcopy unit. C 50 CALL PLTINT C C Determine the axis type C CALL AXTYPE (IQUAD) C C Find the scale for the plot C CALL SCALE (IQUAD,DVSCL, IDVOFF,RIDVSC, IIDVOF) C C Draw the axes. C CALL PAXES ( IQUAD) C C Plot the forces. C DO 400 K = LOWPT, HIGHPT C IX = IFIX ( RIDVSC * RIDV ( (K - LOWPT) + 1 ) ) + IIDVOF IY = IFIX (DVSCL *DV ( (K - LOWPT) + 1 )) + IDVOFF C CALL PLOTER (.IX, IY ) C 400 CONTINUE C CALL PMOVER (0,0) C C Print out the maximum and the minimum values of the force. C WRITE ( IOUT, 200 ) DVMAX, DVMIN 200 FORMAT ('1'/'1'/'1'/ +' The maximum value of the variable is', G11.4, /, +' The minimum value of the variable is', G11.4, C C Pause for a minute to allow the user to observe the maximum and minimum C values of the chosen force. C WRITE ( IOUT, 4100 ) 4100 FORMAT ( /'1' /'1' /'$Press RETURN to continue:' ) READ ( IN, 4200 ) ANSWER 4200 FORMAT ( Al ) 98

Print file "HRDCPY" Page 2 C C Return to the calling program. C RETURN END C C g9

Print file "PRNT" Page 1 SUBROUTINE PRNT C C General Description: C C This routine is used to print out the data from a data file. C The user can choose whether he wants the data sent to the video terminal C or the line printer. C INTEGER POINTS, OPTION LOGICAL* 1 ANSWER C COMMON /INOUT/ IN, IOUT COMMON /ACQBLK/ POINTS,PARAM C DIMENSION PARAM(1000) C CALL ASSIGN (1,' KDX2.DAT' ) CALL AS SIGN (2,' KDY2.DAT' ) CALL ASSIGN(3,' KDZ2.DAT' ) CALL ASSIGN(7,' LP:' ) C C Display a header. C 10 REWIND 1 REWIND 2 REWIND 3 C WRITE (IOUT, 12) 12 FORMAT(' 1'/'1'/'1' / + 29X,'PLOTTING PROGRAM' C C Display the printing menu. C WRITE ( IOUT, 20 ) 20 FORMAT(' 1' / + 33X,'PRINTING MENU' /// + 22X,'1) Send values to the printer' / + 22X,'2) Send values to the screen' / + 22X,'3) Return to the main menu' //////// ) C C Prompt user for menu option. C WRITE ( IOUT, 30 ) 30' FORMAT ('$Please choose option ( 1, 2, 3 )' ) READ ( IN, 32, ERR = 10 ) OPTION 32 FORMAT ( I2 ) C C Check to see that the user chose a valid option. C IF ( OPTION.GT. 3.OR. OPTION.LT. 1 ) GO TO 10 C C Branch to carry out the option requested C GO TO ( 40, 50, 1000 ) OPTION C C Since the user wants to send the force numbers to the printer, C reassign the output device. C 40 IOUT = 7 C C Print a heading on the top of the page. 100

Print file "PRNT" Page 2 C 50 WRITE ( IOUT, 100 ) POINTS 100 FORMAT ('1'/'1'/'1'/ + 17X,'DATA ACQUISTION RESULTS' /// +' Number of samples:', I3,' (unitless)' +' Force values given in Newtons' /// +' SAMPLE #', 10X,'FX', 15X,'FY', 15X,'FZ', /// ) C C Loop to print out the results. C DO 3000 J = 1, POINTS C C Print out the results. C READ (1) M,FX READ (2) M,FY READ (3) M,FZ C WRITE ( IOUT, 2000 ) M, FX, FY, FZ 2000 FORMAT ( 2X, I3, 3 ( 2X, G16.4 ) C C Continue the looping. C 3000 CONTINUE C C If the force values were sent to the screen, pause a moment to allow the C user to view the last of the values printed. C IF ( OPTION.NE. 2 ) GO TO 600 WRITE ( IOUT, 5910 ) 5910 FORMAT ( / / /'$Press RETURN to continue:' ) READ ( IN, 5920 ) ANSWER 5920 FORMAT ( Al ) C C Reassign the output device to be the screen just in case it was changed C to the line printer above. C 600 IOUT = 6 C C Return to the printing menu. C GO TO 10 C C Return to the calling program. C 1000 CLOSE (UNIT=l) CLOSE (UNIT=2) CLOSE (UNIT=3) C RETURN END 101

Print file "CHOICE" Page 1 SUBROUTINE CHOICE C C General Description: C C This routine selects the appropriate variable for plotting. C CALL AS SIGN ( 1,' KDX2.DAT' ) CALL ASS IGN (2,' KDY2.DAT' ) CALL ASSIGN(3,'KDZ2.DAT' ) C INTEGER POINTS, OPTION LOGICAL* 1 ANSWER COMMON /INOUT/ IN, IOUT COMMON /ACQBLK/ POINTS,PARAM C DIMENSION PARAM(1000) C C Display the selection menu. C 10 WRITE ( IOUT, 20 ) 20 FORMAT('1' / + 33X,'SELECTION MENU' /// + 22X,'1) Select FX' / + 22X,'2) Select Fy' / + 22X,'3) Select Fz' //////// ) C C Prompt user for menu option. C WRITE ( IOUT, 30 ) 30 FORMAT ('$Please choose option ( 1, 2, 3 ) READ ( IN, 32, ERR = 10 ) OPTION 32 FORMAT ( I2) C C Check to see that the user chose a valid option. C IF ( OPTION.GT. 3.OR. OPTION.LT. 1 ) GO TO 10 C C Branch to carry out the option requested C GO TO ( 40, 50, 60 ) OPTION C C Since the user wants to send the force numbers to the printer, C reassign the output device. C 40 REWIND 1 DO 21 J = 1, POINTS READ (1) MSTEP, FX D FX = FX + 780.65 IF(MSTEP.GE.17) FX = FX + 21.88 IF(MSTEP.GE.33) FX = FX + 55.88 IF(MSTEP.GE.50) FX = FX + 21.49 IF(MSTEP.GE.63) FX = FX + 13.28 IF(MSTEP.GE.73) FX = FX + 23.05 PARAM(J) = FX 21 CONTINUE GOTO 1000 C 50 REWIND 2 DO 22 J = 1, POINTS 102

Print file "CHOICE" Page 2 READ (2) MSTEP, FY D FY = FY + 217.24 IF(MSTEP.GE.17) FY = FY - 1.96 IF(MSTEP.GE.33) FY = FY + 37.51 IF(MSTEP.GE.49) FY = FY + 17.97 IF(MSTEP.GE.62) FY = FY + 18.37 IF(MSTEP.GE.72) FY = FY + 17.97 C PARAM(J) = FY 22 CONTINUE GOTO 2000 C 60 REWIND 3 DO 23 J = 1, POINTS READ (3) MSTEP, FZ D FZ = FZ + 343.83 IF(MSTEP.GE.17) FZ = FZ - 8.99 IF(MSTEP.GE.33) FZ = FZ + 58.60 IF(MSTEP.GE.49) FZ = FZ + 22.27 IF(MSTEP.GE.62) FZ = FZ + 12.90 IF(MSTEP.GE.72) FZ = FZ + 19.93 C PARAM(J) = FZ 23 CONTINUE GOTO 3000 C 1000 CLOSE (UNIT=1) GOTO 5000 C 2000 CLOSE (UNIT=2) GOTO 5000 C 3000 CLOSE (UNIT=3) GOTO 5000 C 5000 RETURN END 103

Bibliography [1] Lauderbaugh, L. K., and Ulsoy, A. G., "SIMULA - Digital Controller Simulation Package," Technical Report No. UM-MEAM-84-22, Department of Mechanical Emgineering and Applied Mechanics, University of Michigan, Ann Arbor, Michigan. [2] K. Danai and A. G. Ulsoy 1985 in Sensors and Controls for Manufacturing 137-148. A Dynamic State Model for On-Line Tool Wear Estimation in Turning. [3] Takahashi, Y., Rabins, M. J., and Auslander, D. M., Control, Addison-Wesley Publishing Co., Inc., Reading, Mass. 1972. [4] G. C. Goodwin and K. S. Sin 1984 Adaptive Filtering, Prediction, and Control Englwood Cliffs, New Jersey: Prentice-Hall. [5] Rasmussen, F., Lauderbaugh, L. K., and Ulsoy, A. G., "RELTIM - A Library of Routines Used in Real Time Machine Control Applications," Technical Report No. UM-MEAM-84-24, Department of Mechanical Emgineering and Applied Mechanics, University of Michigan, Ann Arbor, Michigan. 104

i9015 02514 786 II THE UNIVERSITY OF MICHIGAN DATE DUE,Ai'z.. /tj'I -} elH Jo, m', 1