<http://www.w3.org/2009/08/skos-reference/skos.html#note>
"A PARAMETER\u2014ESTIMATION ALGORITHM FOR SMALL DIGITAL COMPUTERS ROBERT JAMES TAPP B.SCe, University of Victoria, 1969 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in the Department of Electrical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July, 1972 In presenting this thesis in partial fu1fiIment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of fcc-fncg/ E\"nq in&amp;er in The University of British Columbia Vancouver 8, Canada Date Au.gufr-1 2., /3\"7*X ABSTRACT An algorithm is developed for performing parameter esti mation on a small-size digital computer. First principles of matrix algebra are used to derive a sequential estimator which computes an estimate of a general parameter array A from an array of measurements Z = HA+V where V is a matrix of zero-mean noise terms. At every stage a new row is adjoined to each of Z, H. and V and a new estimate of A is calculated recursive ly, with any one of three well-known filtering processes avail able from the same basic set of recursive equations: aleast-squares filter to minimize J = trace (Z - HA)(Z-HA)-, a maximum-likelihood filter to maximize pZjA(Z|A) or a maximum-A \"\"\"\" *\"\u2014 a-posteriori filter to maximize PA|Z(AIZ). Provision is made for starting the filter either with a-priori means and vari ances of the parameters or with a deterministic \"minimum-norm\" composition based on the first s measurement rows, s being the number of rows in the parameter array. The algorithm is applied to the problem of identifying the parameters of a discrete model for a linear time-invariant control system directly from sequential observations of the inputs and outputs. Results from computer tests are used to demonstrate properties of the algorithm and the important computer programs are included, along with suggestions for further applications. ii TABLE OF CONTENTS Page Xo Int JTOdUC tlOl\"! ooeo\u00A9o\u00ABooooo\u00BBoo\u00BBo*\u00BBo\u00A9o\u00AB\u00AB\u00A9ooQee\u00ABeoooo\u00A9 X II\u00BB Least-Squares Filtering \u00AB . 3 ITIo Maximum-Likelihood Filtering ........ . ..\u201E 20 IV. General Computational Algorithm ....... 29 Vo Identification of A Linear Stationary Process ..oo 37 VX \u00A9 EX!cHTl]p X 6 S oo\u00ABoo\u00ABo\u00ABeooo04oooo*i&gt;oo\u00AB*o**o\u00AB\u00BBooBe\u00ABeooeo\u00AB 43 \u2022Vile 'Further \u2022App-lioartions - \u2022 57 Rapid Identification 8 Identification of Non-Linear Systems o.o..&lt;,\u00AB&gt;. 61 Time-Varying Parameters oo 63 REFERENCES ooft\u00A3&gt;ooooao\u00ABe\u00BB*e*\u00BB*\u00ABcr*\u00ABoooe*\u00ABoe\u00AB4&gt;o\u00ABo*\u00BB\u00BBoo\u00ABe\u00BBft' 65 APPENDIX o*o\u00ABooooooo0eaoo\u00AB\u00AB*fteo\u00ABe*aoo*\u00ABtteo\u00AB0\u00BBooooaoo\u00BBoo 63 Program-Equivalents of Symbols Appearing in Memory-Allocation for Identification Programs 69 Instructions for Using the Identification Pro-Instructions for Writing I/O Subroutines .... 73, Identification Programs for Nova Computer \u00BB\u00BB&lt;&gt; 78 iii LIST OF TABLES Page ...Lo ...Es.sen.tial .D.if.f.er.enc.es ...of ,.L.e.as.t-.Squares., Maximum-Likelihood and Bayesian M.A.Po Filters ..\u00AB.&lt;&gt;...ooo 29 II. Stage-wise Errors of Least-Squares Filter (Example iv LIST OF FIGURES Page ,.l,o ..General .Computational Algorithm .for Estimation 3.0 2. System Identification Algorithm .....&lt;&gt;.... ..00.... 40 3. Time-response of continuous system corresponding to 4o Plot of performance functions for example 2 ...&lt;&gt;. \u201E 49 5. \"Plot of estimation errors for example 3 \"CD '51 6. Plot of estimation errors for example 3 (2) 52 7. Plot of estimation errors for example 3 (3) .0.... 54 8. Plot of errors for estimates in example 4 ........ 56 v ACKNOWLEDGMENT I should like to express my gratitude to the National Research Council for the financial support afforded me in the form of a Post-graduate Scholarship (1969-70) and a Post graduate Bursary (19 70-71), and also to the following indi viduals: my supervisor, Dr. E.V. Bohn, for his patient guid ance; Dr. Ro Donaldson for his helpful suggestions; Dr. Graham Qualtrough for his useful ideas; Mr. Dave Holmes for his as sistance with the computer work; and my good friend and fellow student, Mr. John Wong, for his advice and encouragement. Ro Tapp, July, 1972o vi 1 lo Introduction When it is necessary to estimate important parameters of a system from measurements of system variables, the choice of an optimal mathematical procedure depends on the amount of statistical information available concerning the system and measurement process. Unfortunately, not enough information is available in many practical situations to permit using well-known estimators like the Kalman filter, nor is it obvious how these procedures can be adapted for simpler problems. Kishi [sj? Sage J^l3J, Young [ivj and other authors have indicated how classical least-squares filtering can be useful because of its validity in. the absence of statistical information and -its similarities -with .more sophisticated -me-fcnod-s-,- -but -very little has been written in the way of a unified and complete theory of practical least-squares filtering,, Greville J^J presents a derivation of least-squares curve fitting which is mathematically rigorous but unnecessarily complicated by the use of generalized-inverse theory and not directly applicable to the problem of parameter estimation. In an attempt to apply it to the estimation problem, Kishi J^8J loses some of the mathematical rigour and neglects some important practical considerations. Young [ivj and Sinha and Pille ^15J have con tributed accurate but very simplified descriptions of the method. There is considerable advantage to be gained by using a classical least-squares estimator as the basis for on-line filtering algorithms because it is straightforward to imple-2 merit, valid under most conditions and easily modified for a-priori statistical information. It is the purpose of this thesis to develop a complete theory for least-squares filter ing, leading to an algorithm that can be programmed on a small digital computer and to considerations of how the algorithm can be extended for a number of practical situations\u00AB The mathematical approach used by Greville [2,3 j was chosen as the most suitable on which to base the derivations for general least-squares filtering equations, although his use of Pen rose's pseudo-inverse theory \u00A3ll, 12 j has been abandoned in favour of a more straightforward approach which employs only first principles of matrix algebra. To include the statistical maximum-likelihood and Bayesian filters, some simple modifi cations of the equations are considered. In this thesis all symbols representing vectors and ma trices are underscored, with upper-case letters denoting ma trices and lower-case letters denoting column-vectors wherever possible. A symbol followed by a prime indicates the transpose of the corresponding matrix or column-vector ( example: A )e Where dimensions of a matrix or vector are given, they are enclosed in parentheses following the symbol (example: B(mxn)). M It The identity matrix is represented by the symbol I and matrix tr 11 inverses are denoted by the superscript -1 \u2022\u00AB The symbol for 11 ti the statistical expected-value operator is e . 3 - H\u00BB Least-Squares Filtering An arbitrary but very general representation of the re lation between a collection of measurements of system variables and the basic parameters of the system is HA + V (1) where Z is an array containing all the measured data, A is the array of unknown fixed parameters, H is the matrix representing the defined relationship between the quantities measured and the parameters, and V is an array of measurement noise terms. In a simple example of a body moving with a constant velocity v, it is desired to estimate the velocity and initial position s -of the -bodv -from measurements of its-Dosit-i'on -s -at-known o times t\u201E The parameters Sq and v are defined by the equation s = s + v t o If the position is measured at times t^, and t^ and values s^, S2 and Sg are obtained, then a representation corresponding to equation (1) would be - \u2014 1 fcl *2 =. 1 H *3 1 o V n. n, n. where n^, and are measurement noise terms. The classical method of least squares assumes that for A zero-mean noise the estimate A of the parameter array A should 4 result in a minimum of the sum of the squares of the elements of the matrix (Z - HA), This corresponds to minimizing the cost function J = \u2122 trace (Z-HA) (Z-HA)' (2) If the rows of Z are labelled successively oo\u00AB 9 and \u00AB \u00AB \u00AB the rows of H are similarly labelled \u20142* ~3S*00' tnen tne cost function, can be written J s fp^-h.AKz^^A) (3) For a minimum the derivative with respect to A must be zero: ^ ~ -1~} \" ~ H'Z \u00BB H'HA (4) If the number of rows in H is greater than or equal to the number of columns and the columns are linearly independent then the column vector Hu, which is a linear combination of the columns of H, is non-zero for all non-zero u o Therefore u H Hu is positive for all non-zero u which means that H H is positive definite and hence non-singularo (4) then gives the unique solution 5 A = (H'H)\u21221 H'Z (5) If the number of rows in H is less than or equal to the number of columns and the rows are linearly independent then the row vector u H, which is a linear combination of the rows of H, is non-zero for all non-zero u. Thus u HH u is positive g \u2022f oral! non-zero -u -and -H H is positive -definite- -and therefore nonsingular- Pre-multiplying both sides of (4) by H gives H H'Z = H H*H A A Z = HA (6) Except for the case where H is square, this equation does not have a unique solution, but although no unique solution can be defined on the basis of the least-squares criterion alone it will nevertheless be desirable to define some arbitrary solution,- The most logical choice is that least-squares solu tion which has a minimum \"norm\" and is found by minimizing the cost function J = i trace (A A' ) (7) n &lt;L ------. subject to equation (6)\u00AB Using Lagrange's method of undeter mined multipliers, an augmented cost function is defined: J = trace a \"I A A f A | AA + MZ-HA) (8) where X is the array of undetermined multipliers. Now 6 = A - H X = 0 A a H X (9) Using (9) in (6) Z = H H X X (HH ) (10) Using (10) in (9), A \u00BB H (HH ) Z (11) -Thi-sequation - will - define\"the \"leas t-squarres --e-st-ima-teof A whenever the number of rows of H is less than or equal to its number of columns and the rows are linearly independent.. For the many applications where the observations are not available all at once but are received sequentially in time, it is desirable to have a recursive relation which will provide parameter estimates at every stage by updating prior estimates as each new set or block of data arriveso The addition of more data to the Z matrix will require that elements be added to the H matrix and since the dimensions of H will be changing at every stage it is important to establish which of equations (5) and (11) should be used to determine the estimate at each stage. If the parameter matrix A is to have fixed dimensions, 7 labelled (sxr), then equation (1) shows that H must always have s columns and Z must always have r columns. Thus in this scheme j elements adjoined to the H and Z matrices at sequential stages must take the form of additional rows. If q is the number of rows adjoined to each of 22 and H at every estimation stage, then the total number of rows in each matrix is kq where k is the number of the current estimation stage\u00BB To summarize the dimension labels, (1) can be re-written Z(kqxr) = H(kqxs) A(sxr) + V(Jcqxr) (12) Now, using (5) and (11), the least-squares estimate for A at stage k is defined by A]c = (H^)\"1 Zk , kq&gt;s (14) where Hj, and Z^ are the matrices H and Z at stage k. If qSs then (14) will apply for all values of k, but if q&lt;s then (13) will apply until k exceeds ~ and (14) will apply for all further stages* In designing a general recursive relation for (13) and (14), advantage can be taken of the fact that both solutions would apply for a stage k where kq = s, provided the rows of H are linearly independent. would be square and nonsingular and (13) and (14) would reduce to , A]c - H\"1 Zk , kq = s (15) 8 Thus if the number of rows adjoined to at each stage (q) is a factor of its number of columns (s) then there will be a stage where kq=s such that both (13) and (14) are valid and the final solution from the recursive form of (13) canbeused as the starting value for the recursive form of (14). To obtain the recursive forms for equations (13) and (14) it is convenient to introduce new symbols and defined by Sk = \u00A3k&lt;HkHk) -1 (16) (17) -In the \"theory -of -generaiized invers-es Gk\"would -be \"called the right generalized inverse or right pseudo-inverse of H_k and J_k would be called the left generalized inverse or left pseudo-inverse of The matrices Z^, H_k, G^ and are partitioned as follows: Z^(kqxr) Z^^^ ( [k-l] qxr) 0*000000000000 Z^(qxr) (18) H^fkqxs) = Hkral( ||c-l]qxs) oooooooaoooooo Hk(qxs) (19) G_k(sxkq) = Fk(sx [k-l] q) : E]c(sxq) (20) Equation (13) can now be written .Using \u00A319.)., this can.be written ,as jtwo. .equations.. From (25) Substituting this into (22) gives AK - AK-1 - EKH\u00ABAk,1 + EkZj , kq2s 9 Jk(sxkq) = D. (sx [k-l] q) I B,(sxq) (21) *k - ~k~k = F-Z k~k-l + -k-k ' kq=; (22) ,.T.o .salve .for ,E_k ...and .E_k., define the matrix 2k - 2kHk \u2022- iCciyir1^ - ~k~k-i + h\u00A3 (23) Post-multiplying by Hk gives H, \u00AB F. H. .H. + E, H* H, ~k \u2014k\u2014k-l\u2014k \u2014k\u2014k \u2014k (24) ^k-1 = \u00A3-k\u00BBk~l\u00BBk-l + ^k^k-1 (25) F, H, , Hr + E. Hr H. ~k -k~k~l-k -k-k -k (26) L - Mu_n H..,) A - E, j\u00A3 H,;_, (a ,H,\" ;:) ,-i -k ~k-l ~k~l~k-l k-k -k-l -k-l-k-1' (27) A A \u00A3k = \u00A3)c-l + - \u2014k ~k-l ^ &gt; kc* = s (28) 10 and into (23) gives and into (26) gives Hk =Qk^Iik -\u00A3kH^Qk^ + \u2014k~k \u2014k .Equations...(2.8.),, ..(.29.) and ,(.30.) constitute .the,..rjec..ur,s.i.v:.e..,relation which corresponds to equation (13)\u00BB It may be verified from these equations that the correct starting values for A and Q are zero, for then % = HJ'CHJHJV1 L = H*' (H?Hf V^Z? \u20141 \u20141 \u20141 \u20141 Si 3 H*' &lt;\u00A3* H*' J\"1 which are consistent with the definitions of A. and Q. in (13) and (23)* Using (18) and (21), equation (14) can be written Hk(I&lt; (30) 11 (31) To solve for and , begin by forming the product H^J^. using (17), (19) and (21); R J., = H, (H.'H, J\"1 H' \u2014k\u2014k \u2014k \u2014k--k \u2014k \u2022 -k~i-k eoooooooeoooeoo (32) Pre-multiplying by gives H, S5 H. \u2014k \u2014k H, _ D, \u00BB H, - B, \u2014k-l\u2014k o \u2014k-l\u2014k oooooaoooooooeo ^ I ^k (33) Using (19) 'this can -he \u2022written a-s \"the -two -equatrons (34) Hk - Hk-l-k-l-k + $ (35) Prom (34) 2* \u00AB ^\u2014k-l\u2014k-l + \u2014k' \u2014k^\"\"1 \u2014k-l (36) and from (35) \u2014k = (2k-A-l + (37) If a new matrix is defined by 12 13 D. = P. .H, , - B. H*P. -H. , (44) =k -k-l-k-1 ~k~k ~k-l-k-l and in (40) gives R B p * \u2014 B H* P ' -k -k-l-k -k-k-k-l-k Bk - Pfc^Cl +igPk_1^,r1 (45) Using (44) in (31) Ak - PJc\u201E1H^1Zk_1 - 2k&amp;2k-1l\u00A3\u201E12k-1 + BkZk , kq-s e \u20141 Since \u00A3jc^j_ = (Jl]c_xi-jc_x ^ \u00BB the last equation becomes 4 - + Bk(zJ - HjA^) , kq-s (46) Equations (43), (45) and (46) provide the recursive relation corresponding to equation (14) and can be started by applying (14) and (38) directly to the first stage k such that kqis, which will require inversion of at least an s x s matrix. Since matrix inversion requires fairly complex programming on a small computer, it is perhaps better to arrange that the starting value for (46.) be taken from the last solution of (28) at a stage k where kq ~ s, as described earlier. Similarly a recur sive relation can be found which will provide a starting value for \u00A3k in (43) when kq = s. From (38) the definition of P_k is P, = J. J,* \u2014k \u2014k\u2014k 14 and from equations (16) and (17) 2k = \u00A3k = HjT1 , kq = s Therefore Pk - - GkGk , kq = s (47) Thus at a stage k where kq\u00BBs it is possible to obtain the starting value for from a recursive relation for Rk = \u2014k\u2014k - ^kV\"1 \u00BBk (48) Using (20) this can be written Rk \u00AB \u00A3k!lk + EkEk (49) From (27) \u00A3k = (I \u00BB IkH^H^^H^H^)\"1 (50) Substituting this into (49) gives \u00A3k - (I\"^^^^) \"(H^H^) X (I-EfcHj) + EkEk = (l-^H^R^a-E^)' + EkF^ 15 \u00A3k 51 ^k-i \"* -k-A \u00A3k ~ \u00A3k-i + h^'^l^X + -k-k (51) which is in a convenient form to be calculated in conjunction with (28), (29) and (30). The final general algorithm for least-squares estimation of the parameter matrix A would therefore use equations (28), (29), (30) and (51) for all estimation stages k such that kq 2s and for all subsequent stages would use equations (43), A (45) and (46) beginning with the values of A^ and given by (28) and (51) at a stage k where kq = s\u00AB The calculations involved in these equations are easily performed on a small computer, apart from the following in verses which appear in (30) and (45) respectively: Cl + H^-iH*')-1 As shown in (19) the dimension of is q xs which indicates that both of the matrices being inverted above have dimension qxq, q being the number of rows adjoined to Z and H at each estimation stage&lt;&gt; Thus by choosing q = l, both inverses will involve scalars and the necessary computer programming will be vastly simplified. The number of rows adjoined at each estimation stage need have no effect on the number of rows adjoined at each measurement stage because the measured rows can be stored and adjoined in the estimation algorithm one at 16 a time. Selecting q = 1 also has the advantage that q will always be a divisor of s, the number of columns in H, which is the requirement for proper linking of the two sets of equa tions as previously explained. When q = 1, the matrices and H* degenerate to row vec tors and E^. and B_k degenerate to column vectors. For this reason it is desirable to change the notation \"and replace \u00A3k b* \u00A3k \u00A3k b* ^k Equation (30) now becomes (52) If the column vector (1 - Qk\u201E-j. ^\u2014]c. ^n ^his equation is given the symbol c_k, \u00A3k = (I-WUfc- (53) then from the definition of Qk inequation (23), which was Qk - GkHk . H^CH^)-1^ it can be seen that \u00A3k \u00AB H^I-Qk^)' = hk(I-Qk^) (54) and 17 \u2022 \u00A3k &lt; I - Qk-i - %-i+ ^c-A-i ^k (55) so that equation (52) can now be written (56) and equation (29) now becomes Qk \u00AB Qk^ + ekhk(I-Qk^) ^k-1 + ^k (57) Following is a summary of the major equations and their starting values for the simplified algorithm where q=l: \u00A3k - (I-Qw^ e,_ = c, (c.c,_) (53) kSs &lt; \u00A3k - \u00A3k^k\u00A3k (56) Qk = Qk-! + \u00A3k\u00A3k , 2k = \u00A3 at k = 0 (57) 18 \u00A3k = \u2014k\u20141\u2014k\u2014l~k\u2014k \u2122k\u2014k\u2014k\u20141 + \u2014k\u2014k~k\u20141\u2014k\u2014k + \u2014k\u2014k \u00BB k fs &lt; Rk = 0 at k = 0 A, a 0 at k = 0 (58) (59) k^s &lt;{ \u00A3k - - bkhkPJc_1 , Pk = Rk at k= s -4k = 4k-l +..Sk(^k-^c^l) * (60) (61) A^ ~ A_k from ( 59) at k = s (62) It has already been shown that when the rows of the matrix Hk are linearly independent, the product Hkjik is nonsingular and the matrix is a left identity for the matrix Hk because QkHk - \u2014k * \u2014k\u2014k *\"1 \u2014k\u2014k - % Similarly Qk is a left identity for any other matrix whose columns lie in the transposed row space of Hk\u201E Thus if a new 19 \u00AB row h_k is adjoined to the R matrix and is a linear combination of the previous k-l rows, then the vector h. will lie in the transposed row space of H^_1 and equation (53) will give ck = (I-Q^)^ -= 0 (63) Since the recursive least-squares procedure requires that H have maximum rank at every stage and in particular that the rows be linearly independent for the first s stages, c^, being the first calculation involving a new row of H, is an extremely useful indicator of this conditiono Depending on the process involved, a measurement which would make the rows of H linearly dependent can be rejected in favour of a new measurement or the entire process can be re-started with a minimum of wasted time0 Although at this point all the essential equations for a least-squares filtering algorithm have been developed,a pre liminary comparison with statistical methods will lead to minor improvements which make the algorithm much more useful. In the next chapter will be presented a derivation of the statistical maximum-likelihood filter which parallels that of the least-squares filter in this chapter. Chapter IV will then describe the complete mechanics of the final computational algorithm which was used in the research project outlined in Chapter V. 20 IIIo Maximum-Likelihood Filtering A maximum-likelihood procedure gives the optimum minimum-variance parameter estimate when no a-priori statistical in formation is available concerning the parameters and the noise terms affecting the measurements are zero-mean independent ..white-.Gaussian random ..variables .of .known variance.* The development of the maximum-likelihood filtering equa tions in this chapter follows closely that of the least-squares filter in the previous chapter in order that similarities between the two methods will be apparent. This should facili tate explanation of the general-purpose computational algorithm to be presented in Chapter IV. As in Chapter II5 equation (1) will be the arbitrary representation of the measurement process, where as before the matrix H is assumed to have maximum rank. The optimum estimate A of the parameters A is chosen so that the probabil ity density of each measured quantity conditional on A=A has a maximum at the observed value of the measured quantity. The probability density function involved is often given the name \"likelihood function\" and since the noise terms are statisti cally independent the likelihood function for a row i of mea surements is the product of their individual likelihood func^-tions, which are Gaussian: A (2u) r/2 r/2 exP s \u2022 (64) where r is the number of measurements in a row or the number 21 of columns in the measurement array and it has been assumed that the noise on each measurement of a row has the same vari ance s^o This latter assumption results in a great simplifi cation to the derivation which follows and does not seriously limit the usefulness of the equations, because measurements having different noise-variances can always be located in separate rows. A product of the likelihood functions of all k rows gives the likelihood function for the entire measurement set at stage k: P2\u00ABA=(27r)kr/2(de77~r72exP IT\"*\u2014 1 ' * ' ' \u00BB \u2022g-\u00ABL s^' (- hu A) (- h-A) (65) where S = sx 0 0 0 s2 0 0 0 s0 , a positive-definite matrix* Maximizing this likelihood function is equivalent to maximizing its logarithm? 1T3 \u20141 , \u00BB \u00BB \u00BB t \u00AB tr log pz{A = -2^-si ~ 111 A) (\u2122i ~ ~i~) ~ ~T log (27t ) - | log (det S) (66) and a maximum results when the derivative with respect to A is zero: 22 H'S\"1 Z = H'S\"\"1 HA (67) When H has fewer rows than columns, (HH ) exists and the last equation reduces to Z = HA (68) This is identical to the least-squares result of equation (6) and the minimum-norm estimates for the two methods are the same: A e H* (H H' J\"*1 Z (69) \u00BB When H has more rows than columns, H H is positive definite 1 \u2014 1 and so is H S Ho Thus the maximum-likelihood estimate from equation (67) is A = (H'S\"\"1 H)\"1 H'S\u21221 Z (70) which is the least-squares solution of (5) weighted by the inverse of the noise variance matrix S_s A recursive relation for equation (70), unlike the method of least squares, cannot theoretically be started using the minimum-norm result of (69) at stage k = s because (69) does not contain the information regarding the noise variance for stages k \u2122 s that is required by (70). This problem will be discussed later. It is first necessary to obtain a recursive 23 form for (70K Following a procedure similar to that used for the re cursive form of (5), the maximum-likelihood estimate at a stage k is written as (71) where J_k is now defined by (72) To make the equations compatible with the algorithm derived near the end of Chapter II, one row only will be adjoined to each of and at every stage, and the following partition-ings are valid: Zk(kx r) t 4c ooooooeooooeo z, (lxr) Hk(kx s) = Jk(sx k) -k-l( fc\"\"1] x s} I hk(1x s) Vsx[k-l] &gt; I bk(sxl)j (73) (74) (75) (j \u20141 O -k-l I 0 oeo\u00AB\u00ABoooooe -1 (76) Now 24 oeeo\u00ABooooeooAOO h'p. ^A (77) ' \u20141 Pre-multiplying by H^S^% \u2022 -1 H, S. \u2014k-~k ' -1 \u00A3A -k-i-k : %-A HA : hkbk (78) Using (74) and (76) this can be written as the two equations Sk-A-i - \u00BBk-A-A-A + Vk1 HA (79) k1 - \"k-A^A-A + Vk\"1 HA (80) From (79) \u2014k A-A-A-i + Vk1^\"1^!^! (81) From (80) *k - &lt;\u00BBk-A-A-i + Vk1 V\"1 Vk1 (82) Defining \u00A3k - AS-1^)-1 = (H^S-^H^ + h^1^)\"1 (83) (81) and (82) become -k ~k-k\"k From (83) Pre-multiplying by P. and post-multiplying by P, Using (85) this becomes and using this result in (84) gives and in (85), 25 2k = \u00A3kHk~Ai (84) b, = P. h, s, 1 (85) ^k1 \" \u00A3k-\\ + Vk\"1^ (86) *k-l = *k + ScVk * 4-1 (87) *k \u00BB \u00A3k-l - ^AA-1 (88&gt; 2k \" Sk-X-iC1! \" ^4-iCi4i (89} kk - Ek-iVk1 - fek^A-iVk'1 bk - \u00A3k-A(sk + HA-iV\"\"1 (90) Using (73), (75) and (89), equation (77) can be written k - 4-1 + ^k^k-Hkik-^ (91) 26 Comparing equations (90), (88) and (91) with (60), (61) and (62) respectively, it is seen that the maximum-likelihood ing equations except that the \"1\" in equation (60) has been replaced by the variance term s^ in equation (90). In other words the maximum-likelihood filter degenerates to a least-squares filter if s, = 1. The starting values for and can be obtained by a direct application of (70) and (83) to the first s measurement stages, which would require inversion of an s xs matrix* How ever, since starting values constitute a-priori known statis tics of the parameters it is instructive instead to compare the recursive maximum-likelihood filter with a similar filter that is based on such statistics. In the maximum-a-posteriori (MAP) filter, A has a normal or Gaussian probability distrib-A ution, A^ is its expected value and is a diagonal matrix such that the ith element on its main diagonal is the variance of every parameter in the ith row of A. To obtain this filter it is necessary to maximize the so-called a-posteriori density PA|Z which is related to the likelihood density PZ\\A B^ THE Bayes rule: This is equivalent to finding a maximum of its logarithm: filtering equations are identical to the least-squares filter-(92) log p AlZ = log p ZIA + log pA - log pz (93) 27 Differentiating with respect to A results in HA\" PAIZ = SV1 &lt;z-HA) + ^ log PA (94) where log pz is aero because pz is not a function of Ae The a-priori density pA can be written PA = (27t)rs/2(det PQ)r/2 y P_I (A, . - A ^% xi i x, j (95) and then -2- log pA \u00BB -P~ (A-A ) dA 3 ^A o \u2014 \u2014o (96) The maximum a-posteriori density occurs when -g^ log PAj?3 is zero: H'S\"1 tZ - HA) P (A \u2014 A ) = 0 \u2014O ~ \u2014O \u2014 A = P(H S Z + P^ A^) (97) where \u00BB \u2014i \u2014i \u2014i P s (H S H + P^ ) (98) Comparison with equation (86) shows that the MAP estimate will in fact be generated by the recursive maximum-likelihood equa-tions when A and P are the expected value and variance, \u2014o \u2014o respectively, of the parameters. The resulting filter is ac tually a special case of the well-known discrete Kalman filter 28 but. has limited application possibilities because of the need for accurate a-priori statistics of the parameters and because of the restriction that parameters in the same row of A must have the same variance.. However, the fact that the maximum-likelihood estimate of (70) and (83) is generated by the same recursive relations as the MAP estimate of (97) and (98) and that the MAP estimate degenerates to the maximum-likelihood estimate as P becomes infinite, indicates that the maximum-\u2014o ' likelihood filter can be started with A equal to zero and P \u2014o \u2014o a diagonal matrix with large elements on the main diagonalo In this way the recursive maximum-likelihood filter would presuppose A to have zero expected value and very large vari ance, which is consistent with a total lack of a-priori sta-tisticso 29 IV. General Computational Algorithm It is now apparent that with very minor alterations the basic least-squares recursive equations of Chapter II (53-62) can perform'either maximum-likelihood or Bayesian (maximum-a-posteriori) filtering. By substituting the noise variance in .place of the \"1\"- in equation (60) and replacing the minimum-A norm equations (53-59) with initial values A^ zero and very large, a maximum-likelihood filter resultso A Bayesian filter is produced by using the expected value and variance of A for A A^ and respectively in the maximum-likelihood filter. The following table summarizes the differences: TABLE I - -Essential Differences of Least-Squarres, -Maximum- Likelihood and Bayesian MAP Filters Filter s^ Initial Values T _ r. i minimum-norm composition Least-squares 1 . , . Ir-^ m\\ ^ usxng equations (53-59) A Maximum-likelihood noise variance zero, very large A A = expected value of A Bayesian MAP noise variance Zr . _ \u201E x P = variance of A \u2014o \u2014 On the next page is presented a general computational algorithm which allows for any of the combinations in the above table. It also allows for unclassified combinations such as one in which the noise variance term is \"1\" and the A initial values are A zero and P large. This is effectively T\"0 \"\"\"O 30 FIGURE 1 - General Computational Algorithm for Estimation BEGIN ^ YES GET STARTING VALUES P , A CALCULATE bfc, \u00A3k  FROM (60), (61) GET MEASUREMENT z, AND VARIANCE s CALCULATE ESTIMATE A Ak FROM (62) T k = k + 1 (END YES &gt;kMAX^N0 Q3c, Sk--. \u00A3k, Ak - 0 J = S \u2014iz: DEFINE MEASUREMENT BY h, \u2014k CALCULATE \u00A3k FROM ( 53 ) YES CALCULATE e, , Q, , R, \u2014k' \u2014k* \u2014k FROM (56), (57), (58) GET MEASUREMENT z. CALCULATE ESTIMATE A FROM (59) k = k + 1 5 = 5-1 YES 31 a least-squares filter which is begun in the same way as the maximum-likelihood filter, eliminating the need for the \"min imum-norm\" equations (53) and (56-59). Tests of this filter are described in example 1 of Chapter VI. In the general algorithm, linear dependence of the rows of H in the first s stages of the least-squares filter results in a value of c^ which is near zero, re-initializing the entire process. How close to zero c^ must come in order for this to occur is a difficult matter to define and depends among other things on the precision of the calculations. While exact linear dependence would theoretically make c^ exactly zero, the value determined by the computer will normally contain errors due to truncation and thus be slightly different from zero. In any event, cases of near linear dependence can produce inaccurate estimates, so it is probably best to require that c^. remain reasonably large. This can be done by defining a threshold value and causing re-initialization if the magnitude of c, falls less than this threshold during the first s stages. In choosing which of the various filtering procedures to use it is important to know how the errors of the estimates are expected to compare. A useful matrix which gives an esti mate of the error is the error covariance matrix, hereby de fined as A A A | cov (A) = e(A - A)(A - A) (99) 32 The trace of this matrix is equivalent to the expected value A of the sum of the squares of the error matrix (A - A)0 Using equations (5) and (1) it can be seen that the error covariance of the least-squares filter after s stages is given by cov -(ALS-) \u2022= -(-H--H-)\"1 \u2022\u2022H'-\u00A9(V-V^)--H-(-H-'\u00AB\u2022)\u2022\"1 (100) E(VV ) can readily be shovm to be a diagonal matrix such that the ith element on its main diagonal is the sum of the vari ances of the measurements in the ith row of Zo The least-squares estimate is always unbiasedo That is, A e(A-ALS) =o (101) The maximum-likelihood estimate (70) has an error covari-ance given by cov (AM )' = (H'S\"\"1 H)\"1 H'S\"1 e(v v') s\"1 H(H'S-1 H)\"\"1 (102) Because the maximum-likelihood filter requires that measure ments in the same row of zz have the same variance and since there are r measurements in each row, it is apparent that e(V v') = r x S (103) where S_ is the measurement-noise variance matrix as defined in (65)c Therefore, when all measurements in the same row of 33 Z have the same variance, the error covariances of the least-squares and maximum-likelihood estimates become cov (ALS) = r x (H'H)\"1 H'SH(H'H)\"1 (104) cov (AML) = r x (H'S-1 H) (105) The definition of the P-matrix for the maximum-likelihood filter as given in equation (83) shows that the error covar-iance of the maximum-likelihood filter is given simply by cov (AML) = r x P (106) Using the matrix inequality M'M (N'M) ' (N 'N)\"1 (N'M) (107) (see Sage and Melsa [l4J, pa 246)-where M and N are any two kx s matrices with k|s and N of rank s, and making the sub stitutions M. = S 2 H(H'H)\"\"1 (108) N = S 2 H (109) it is easily shown that cov (A_ML) If cov (A ) (110) 34 That is, the maximum-likelihood filter? when applicable, gives, an estimate which is as good as or better than that of the least-squares filter., Like the least-squares estimate, the maximum-likelihood estimate is unbiased: e(AnL-A) .== 0 (111) The error covariance of the Bayesian MAP estimate, as defined by equation (97), is cov (A MAp) = \u00A3 (H'S\u21221 E(V/ ) S\"1 H + Pj\"1 G(AA') Pj\"1) L (112) Since all measurements in the same row of Z must have the same noise variance and the probability distributions of all pa rameters in the same row of A must have the same variance, E(W') a r x S (113) G(A A' ) = r x P (114) where r is the number of elements in each row of Z and A The error covariance therefore becomes cov (A^p).^ r x P (H'S\"1 H + Pj\"1) P = r x P (115) which can be readily determined from the P-matrix. Comparing the values of P for the maximum-likelihood and MAP estimates, it is obvious that the MAP estimate, where applicable, has an 35 error covariance which is less than or equal to that of either the maximum-likelihood or least-squares estimates. In fact, the MAP estimate, when valid, is known to have the least error covariance of any known estimate. Even when the restriction is removed that the noise and parameters be Gaussian the MAP filter still provides the best estimate of all linear filters\u00A9 The noise must still be random with zero-mean and known vari ance and the expected value and variance of the parameters must still be known. The filter is then usually called a lin ear-minimum-variance filter. In addition to the fact that all parameters in the same row of A must have the same variance, the MAP estimate has another major disadvantage. If incorrect prior expected values and variances are used the estimates will be biased, with the bias at a stage k given by e(ik-^} - V^^H^A) + p;1^ - p-^u)) = &lt;I + \u00A3o J^i*!1^\"1 i-e(A)) (116) The bias is most noticeable for smaller values of k and de creases as k increases. It is also smaller for higher values of the initial variance P and approaches zero as P becomes \u2014o \u2014o 36 infinite, the estimate then becoming a maximum-likelihood estimate. In conclusion it may be said that among the three filters, least-squares, maximum-likelihood and Bayesian MAP, the more extensive the a-priori statistical knowledge of the parameters and measurement noise, the lower is the covariance of estima tion error. 37 V. Identification of A Linear Stationary Process The computational algorithm of the previous chapter can be used to estimate the parameters of a discrete model for a linear time-invariant process. If measurements of the system variables are available at uniformly-spaced intervals of time, ,it ..is ...possible to ..develop a .model of the .form where Xjc an n-dimensional vector composed of the system outputs at stage k, u^ is an m-dimensional vector composed of the system inputs at stage k and Q and A are matrices composed of the constant parameters describing the process, x^. is called the state of the system at stage k, u^. the control and &lt;\u00A3_ and A the state-transition and state-driving matrices respectively. Transposing both sides of the last equation results in \u00A3k 3 -k-1^' + ~k-l(118) Because of measurement noise, there will be differences between the observed values of the variables and their true values. Therefore it is convenient to distinguish the observed values with a superscribed bar as follows: Hk = xk + (119) Hk - Hk + &lt;\u00A3k (120) 38 where and are vectors comprised of the noise terms. Combining these two equations with (118) yields the relation between the parameters and the observed values: or x. \u2022k -k-i t -k\u00BBi o e \u00A9 o A + he ~ iik-1^' - 2k-i^ (121) If the measured vectors x^, k = 1,2,3,\u00AB\u00AB&gt;. become the successive rows of the Z matrix in the computational algorithm: z-k (m) (122) and the corresponding prior measurements become the successive rows of the H matrix: iik - \u00A3k-i : Hk-i (n + m) (123) then in accordance with the representation of equation (1) the unknown parameter will be A = o o o o f A (124) and the sequential noise vectors will be 39 xk = u'k - Hk-i^ \" ^k-i^' (125) In other words, use of the z^. and h^ vectors as defined by (122) and (123) in the general computational algorithm of Chapter IV will produce an estimate of the matrix defined by (124). If the components of the noise sequences and have zero means and are Gaussian, white and independent, then the least-squares filter of Chapter II applies because the overall observation noise v^. defined by (125) has zero-mean components. The maximum-likelihood and Bayesian filters, however, are not strictly valid as they have been derived in Chapter III, be cause the components of v^ are not likely to be independent or white. None the less it would seem logical that the hier archy among the three filters should still exist because of the differing degrees of a-priori information utilized. Thus, although methods exist by which the maximum-likelihood and Bayesian filters can be made optimal (see, for example, Sage and Melsa ^14 J, Chapter 8), they involve such extensive com plication of the algorithm that it is convenient in this ap plication to merely ignore the fact that the noise components may be non-white or statistically dependent. The computational algorithm as applied to the system-identification problem is represented in the flow-chart on the following page. All of the experimental tests described in the next chapter were made using this algorithm on either the 40 FIGURE 2 - System Identification Algorithm ^ BEGIN ^ GET SYSTEM DIMENSIONS \u2022GET STARTING VALUES i =1 | YES^STARTINC VALUES 1 = 0 USER'S PRELIMINARY MEASUREMENT SUBROUTINE USER'S SUBROUTINE TO [GET STATE MEASUREMENT k = 1 CALCULATE bk? Pk USER'S SUBROUTINE TO GET STATE MEASUREMENT SUBROUTINE TO CALCULATE ESTIMATE USER'S SUBROUTINE TO OUTPUT ESTIMATE k = k + 1 A J \u00BB s \u2014i\u2014 A,k? \u00A3ks Qkj Sk = o CALCULATE c. = (I \u2014 Q, . ) h, \u2014k , \u2014 \u2014K-l \u2014K CALCULATE e_k, Qk, Rk USER'S SUBROUTINE TO GET STATE MEASUREMENT SUBROUTINE TO CALCULATE ESTIMATE USER'S SUBROUTINE TO OUTPUT ESTIMATE 41 I.BcM, 360-67 or the Data General Corporation \"Nova\" digital computer. In the appendix are included the complete programs for the Nova version of the algorithm. Comparison with the flow-chart of Chapter IV will show that the identification algorithm is basically the same except for the addition of certain specialized subroutines for handling the input and output data. These data subroutines can be changed to suit any particular application. There is a preliminary measurement subroutine which is provided in case there are any tasks as sociated with the measurement process which must be performed before entering the identification cycle. For example, should it be necessary to take samples of the system at a rate faster than the computation cycle v/ould allow, the measurements may all be made in advance and stored by this subroutine. Then on each cycle is a subroutine to get the state measurement from the appropriate source and another to output the calculated estimate*. The algorithm was programmed on the system-360 to allow for more sophisticated analysis using artificial models of known statistics. In the Nova programs, all the initializing procedures are controlled by the operator using the teletype keyboard in a conversational manner. The system dimensions can be set to estimate any matrix A up to a dimension of 8 x 8 and the cal culations are performed in floating-point arithmetic that is based on a 24-bit mantissa with sign and 7-bit exponent using the standard basic floating-point software provided with the 42 computer. Although in the flow-chart the k-counters are separate from the subroutines, in the Nova programs the job of counting stages has been left to the user-supplied subroutines. This allows for counting either at the point where data comes in or at the point of outputting the estimate, whichever is more suitable to the particular application. Also left to the user subroutines is the task of determining the sampling times, which, in the case where data is obtained from an actual system using an analogue-to-digital converter, could require an ex ternal real-time clock connected to the input/output bus of the computer. The programs are thus very versatile, with the permanent software performing identification only, and the user-supplied subroutines having full control over the rest of the process. 43 VI. Examples All the examples described in this chapter were used to test the computer programs for the system-identification al gorithm and to study various properties of the algorithm. Sequential values of the state were generated for the algorithm .in -each case -using. .known values ..of and A .and -the estimates were then compared with these known values. Example 1: The model parameters were 0.995 0.5 0.0 0.0 1.0 0.5 0.0 -1.13 0.9 A = 0.0 0.0 1.25 with initial state and control x = \u2014o 0.0 1.5 3.95 u = 1 The control was left constant throughout the process, resulting in an open-loop response corresponding to the three curves of Figure 3. The curves are shown to be continuous because it has been assumed that in practice the measurements would result from uniform sampling of this continuous system. The simulation was performed using the Nova programs and there was no measurement noise added to the model. Table II 0 12 3 4 5 6 7 8 9 10 11 12 13 M 15 16 17 Figure 3 - Time response of continuous system corresponding to Example 1. TABLE II - Stage-wise Errors of Least-Squares Filter (Example 1) Figures in each column represent errors at successive stages. Stage \"minimum norm\" \u00A3o \u2022 10 x 1 = 102 x I P - 103 x I \u2014o \u2014 P \u00AB 105 x I \u2014o \u2014 P = 107 x I \u2014o \u2014 P - 1016 X I \u2014o \u2014 1 + .4956009E+01 + .4956044E+01 + .4956009E+01 + .4956009E+01 + .495601 1E+01 + .49560S9E+01 + . 4956009E+01 2 +.2193356E+31 + .2197610E+31 +.219340QE+0! + \u20222193356E+01 + .2193354E+01 +.2193356E+01 +.2193356E+01 3 +.7292265E+03 + . 1245663E+31 + .747 6390E+00 + \u2022 729 4487E+00 + .7292240E+00 + .72922 60E+00 +.7292259E+00 4 + .331 5S34E-C3 +.3763182E+00 +.72S6461E+00 + \u2022 2720355E+00 +.2375113E-03 . +.6260436E-05 +.8251297E-0S 5 +.9723237E-03 +.3960244E+03 +.4370673E+00 + .2733907E-01 +.S037871E-05 +.3311357E-06 + .2695410E-06 6 + .H47506E-07 +.S040013E+00 + .1357909E+00 + \u20222373371E-02 . +.S270537E-06 + .6163S93E-07 +.1110659E-07 7 \u2022.1032932E-03 +.S403006E+00 +.3114320E-01 \u20224229341E-03 + .1521597E-06 + .1 694626E-07 + .10322S5E-07 8 +.6331324E-03 +.2660201E+00 +.7173627E-02 + \u20228230054E-04 + .\u2022009 6006E-07 \u2022 .243 IS 67E-0S +.303 6397E-07 9 +.117 42 60E-07 + . 1069380E+00 + .18 60306E-02 + . 1993976E-04 + . 1 173606E-07 +.1211445E-09 \u2022 +.2235941E-03 10 +.1237370E-07 +.4123137E-01 + .569 53 62E-03 + \u20225937541E-05 +.3443537E-03 +.3527151E-09 +.2407231E-B3 11 + .1 1 13249E-07 + .1753404E-01 + .2159944E-03 ; + \u20222214477E-05 + .501435SE-09 + .6 3 7 8 63 0E-G9 +.4007216E-08 12 + .3576017E-0S + .914S61SE-02 + .1058297E-03 + \u2022 1067259E-05 + .2560367E-10 +.1017675E-03 + .2451573E-08 13 + .5319554E-03 + .5936951E-02 +.6645077E-04 .6555694E-06 + .9559106E-09 + . 1640657E-03 + .201S661E-G3 14 + .4583361E-G3 ' +-4574932E-02 +.5036433E-04 + .49 2 AS 7 IE-06 + .3036121E-03 + .22507GSE-03 +.1039153E-03 15 + .4319357E-03 +.3387912E-02 + . 4249 733E-04 + \u2022 419S029E-06 + .1093 628E-07 +.1374427E-07 + . 1715590E-09 16 \u2022 . 432 47 63 E-03 +.3433026E-S2 +.3737925E-04 .37S7415E-06 +.1721376E-09 +.2004652E-03 + .4943047E-03 17 +.4236413E-08 +.30337S9E-02 + .3293429E-04 + \u20223446333E-06 + .3334956E-03 +.2594733E-03 +.2070682E-06 13 + .4066631E-03 +.2629245E-02 +.2355703E-04 . + \u20223090746E-06 + .7399028E-0S + .3 7 2 68 39E-03 +.5057110E-G3 19 + .3959733E-03 +.2224612E-02 + .2403305E-04 + .2663924E-06 +.3632550E-07 +.1630333E-07 +.1006140E-07 on 46 shows the actual estimation error as defined by Error (A^) = trace (A - A^) (A - A^.) (126) computed on the Nova at sequential sampling times for the identification algorithm when used as a least-squares filter ..with .a \"minimum .norm\" ..comp.o.siti.on ..at ..the .start and ...also when started with A^ = 0 and equal to various scalar multiples of the identity matrix 1\u00B0 It can be seen that when is fairly large, estimates can be obtained which are as good as, and at some stages marginally better than those obtained when the \"minimum norm\" procedure is used. Here the filtering problem is a deterministic one because there is no a priori information and the estimates should\"be based solely on the noise-free observations. P^ must therefore be made large to give minimal weighting of the initial estimates A^. The resulting filter is then a good approximation to the purely deterministic least-squares filter of Chapter II, with much less computational requirements. However, the pure least-squares filter is subject to minimum initial bias and with it a better estimate results after fewer measurements* Specifically, the estimation error at the fourth stage in this example was lowest with the pure least-squares filter, the estimates being .9949869 .5000086 -.0000036 .0000224 .9999967 .5000123 -.0000045 -1.130009 .8999898 ^4 .0000246 -.0000396 1.249999 47 Example 2: The model was 0o995 0. 5 0.0 0.0 &lt;L = 0.0 1.0 0.5 A = 0.0 0.0 -1.13 -0.9 1.25 with initial state and initial control ~~o 0.0 1.5 -1.45 u = 1 o and no simulated measurement noise. In generating the remaining states, the control was left equal to Uq until stage 1, after which each state was deter mined using a control chosen to minimize the estimated simple performance function Jk+1 = W+lEi^c+i + HkH2uk where and are weighting matrices chosen for stability A purposes and xk+i defined by the equation k+1 k -k In this example it is possible to have 48 Setting the derivative of Jk+1 with respect to equal to zero gives A j A ^k ^k+1 + Uk = 0 Figure 4 shows the variation of the resulting performance function J(t) = X, X, + U, T , t, ft&lt;t, , \u2014k ~k k~l ' k \u2014 k+1 for the open-loop case where the control was left equal to UQ for all stages and ,for two cases where .u^.was .calculated., beginning with u^, based on estimates' from the least-squares filter using different starting \"procedures. All calculations were performed on the Nova. This method might be useful for simple combined identi fication and control of an actual continuous system by calcu lating a sub-optimal control based on the discrete model. 4. _ J OPEN-LOOP CASE ' FOR \u2022 ALL k - _ J r CONTROL BASED ON SUB-OPTIMAL LEAST-SQUARES IDENTIFICATION J STARTING WITH \u00A3 = 0 AND P = 10xl CONTROL BASED ON LEAST-SQUARES IDENTIFICATION STARTING WITH EQUATIONS (53) AND (56 -59) ) 2 4 8 10 12 14 15 Figure 4 - Plot of performance functions for example 2. J(t) = 4 ^ t uk-l' &lt; t &lt; t k+1 50 Example 3: The model was OoO OoO 0.9 0.0 fL - 2.0 OoO 0.0 A = O.O OoO 0o7 0.0 0.0 with initial state and control x = \u2014o 2.7 10.0 4.9 u = 0.0 The simulation was performed on the system-360 with random -noise -of normal -distri-bution --added\u2014to-\"the-state and -control measurements. The standard deviation of the noise was 0. 7, the variance 0.49. The identification algorithm was used as a \"best case\" A A of the Bayesian I4AP filter, with ^ = j# and A^ = A. s^ at every stage was set equal to the noise variance, 0.49. While in example 1 the initial estimate was inaccurate and the measure ments were exact, in this example the initial estimates are exact and the measurements are noisy. Figure 5 shows the com puted estimation errors as defined by (126) at each stage. As expected, the results are opposite to. those of example 1, with a lower P now giving the better estimates because of increased weighting of A^. The results of Figure 6 were obtained with this same example and show the effect on the estimation error of using 52 53 different multiples of the noise variance for s^ in the al gorithm. When P is large the effect is not noticeable but when PQ is small, increasingly higher multiples give increas ingly better estimates in the stages following stage 4. A higher value of s^. provides decreased weighting of the noisy measurements and increased weighting of the good initial es timate, s^. has no noticeable effect on the estimates prior to stage 4. Figure 7 was obtained by the same procedure, except that the minimum-norm composition was used at the start. 'As-is:,the case when the filter is started with P_Q large, there was neg ligible difference of the errors when values of s^ ranging from 0.5 to 2 times the noise variance were used. Figure 7 - Plot of estimation errors for example 3. Minimum\u2014norm start. 55 Example 4: The model was 0.99 5 0.5 0.0 0.0 0.0 1.0 0. 5 A = 0.0 0.0 -1.13 0.9 1.25 with initial state and control x = \u2014o 0.0 1.5 3.95 u = 1 (see Figure 3 for response curves). The simulation was done on the system-360, -introducing Gaus-sian 'noi\"se of standard deviation 0.5 and variance 0.25 to the measurements. s, in the algorithm was set equal to the variance at each stage. Figure 8 shows the computed estimation errors as defined by (126) at each stage for three different starting procedures: A A Mo=^&gt; ^-^.i \u00A3o=i for a \"best-case\" Bayesian MAP filter; a \"minimum-norm\" composition for a least-squares or maximum-A A g likelihood filter; jZ^ = 0, A^ = 0, P^ = 10 xl for an approx imate least-squares or maximum-likelihood filter. The results still support the hierarchy of filters de veloped in Chapter IV despite the fact.that the overall noise terms defined by (125) are not expected to be statistically independent or white as discussed in Chapter V. 56 Figure 8 - Plot of errors for estimates in example 4. 57 VII. Further Applications What has been presented in this thesis is an algorithm to estimate the parameter matrix A in the general measurement process of equation (1): Z = HA + V (1) While the accompanying computer programs (see Appendix) have been written for the particular system-identification problem of Chapter V, it is a simple matter to adapt them for any measurement process defined by (1). Specifically, the identi-fication problem of Chapter V requires that each row of Z (z^) should be taken from the state measurement which will comprise the next row of H\"('h \u2022[\u2022)\u00BB and 'thus \"for the computer programs the vectors z and h can share the same storage locations. In other applications, separate sets of storage locations may be required. Apart from this, the program, when supplied sequen tially (via the user's measurement subroutine) with the rows of Z and H arrays satisfying the relation of equation (1), will generate a sequential estimate of the parameter array A, subject to the following conditions developed in the previous chapters: The elements of V must have zero expected values. If all elements of the same row of V (v^) have the same probability distribution and the variance of their dis tribution is known, it should be supplied for the value of s, corresponding to that row (maximum-likelihood 58 filter). If all elements of the same row of V do not have the same probability distribution or if the variances of the distributions are not known, Sy should be set equal to 1 at every stage (least-squares filter). If expected values for the parameters are known then A they should be used as the elements of A .If in addition J. \u2014o the probability distributions of all parameters in the same row of A have the same variance and the variances of the distributions of all the parameters are known, then P should be a diagonal matrix with the ith element \u2014o 3 on its main diagonal equal to the variance of the distri bution of every parameter in the ith row of Ao Otherwise P should be a diagonal matrix with each element on its main diagonal set large enough to allow for any uncer-tainty in the corresponding row of A^. If expected values for the parameters are not known A then no initial values should be supplied for A^ and P^, but equations (53) and (56 - 59) should be used at the start. However, where the increased computational time required by these equations would be prohibitive, a good approximation can be achieved by using initial values A A =0 and P diagonal with larqe elements on the main -r-O \u2014 \u2014O 3 diagonal. Rapid Identification A major difficulty with the method of system-identifica tion developed in Chapter V is that the estimated discrete 59 model cannot be accurate unless the rate of sampling the state is high in relation to the rate at which the state varieso At the same time, such rapid sampling can lead to near linear-dependence in the state measurements and consequent ill-con ditioning of the H matrix, which makes adequate identification impossible. Hanafy and Bohn |^4J have suggested augmenting the state measurement at each sampling time with the measured outputs of integrators cascaded to the inputs and outputs of the continuous system. It is claimed that this additional data is effective in overcoming the problem of ill-conditioning. However, the usual treatment becomes cumbersome because the data and parameters must be structured into lengthy vectors in order to fit the form of conventional estimators, which are derived for a measurement process of the type z_ = Ha + v z being the data vector, a the parameter vector, H the mea surement matrix and v a vector of noise terms. For the iden tification problem this results in a large H matrix of block-diagonal form and containing many zeros. The algorithm of this thesis can be used quite readily for identification with augmented state measurements. At each sampling time, the inputs and outputs of the system are mea sured and stored, along with the outputs of the successive integrators. A state measurement is processed as one row in the estimation algorithm, followed by the integrator outputs as subsequent rows. Suppose, for example, that each of the 60 system inputs and outputs is passed through two integrators. Evidently the integrals t _ (t) dt and 0 0 (t) dt dt will satisfy the same linear differential equation as does x (t), so uniform samples of their outputs should satisfy the same difference equation: 2\u00A3 (t]c) A ! x (t) dt = 0 t k-l x'(t)dt k-l u*(t)dt A tr. r t. 0 x (t) dt dt = I 0 I, x (t) dt 0 dt ' \"k-l u (t) dt dt A The beginning rows of data for the estimation algorithm would therefore be t ? \u20141 = \u2014 (t^) X (tQ) a U (tQ) 61 where the superscribed bars indicate that these are the ob served values of the variables concerned. With this procedure the state measurement defining z_ at a given stage does not immediately become the h vector for the following stage as it does in the simple identification problem. Therefore separate sets of memory locations are needed for the z and h vectors, as mentioned at the beginning of this chapter. Identification of Non-Linear Systems Both the identification methods discussed thus far have assumed a linear model for the system being measured. However, it is equally possible, within the allowable forms of measure ment processes, to assume certain non-linear models. Netravali and de Figueiredo J^9J have discussed methods of obtaining re gression functions for classes of discrete non-linear systems in which the evolution operators can be represented by algebraic or trigonometric polynomials. Although noise considerations are more involved, the computational requirements are not unlike those of the linear identification problem and are adaptable to the computer programs contained in this thesis. As a very simple example, suppose it is desired to estimate a third-order non-linear algebraic model of the form Xk+1 = V+ alXk + a2Xk + a3Xk from measurements of the scalar variable x^, k =0,1,2,... . The estimation algorithm would begin with the following data: where again the superscribed bar is used to denote measured values. It is useful to assume non-linear models in some cases involving linear systems where not all of the state variables are measured. For example, although Figure 3 in the previous chapter describes a linear system of 3 outputs and 1 input, a model for any one of the outputs, obtained from measurements of that output alone, would have to be non-linear. Of course, non-linear models are not always necessary to reduce the order of a linear system, because many linear systems can be realized in terms of reduced linear models. A further sophistication of the identification algorithm for linear systems could pro-63 vide for appropriate selection of the measured variables to effect such a reduction. Time-Varying Parameters Time-varying parameters can be accomodated by modifying the algorithm so that prior estimates are updated at each \u2022stage to allow -for expected \u2022time-variations during \"the mea surement interval. That is, if the parameter array A is known to vary according to the difference equation -k+1 = \u00A3 (k+1\u00BB k) ^k then in the algorithm is replaced by its a priori update: = ^k-l This is a much-used procedure and forms the basis of the Kalman filter for state estimation. Other methods are available if no model for the parameter variations is known (see, for ex ample, Young [l7] ). 65 REFERENCES lj Flores, I., The Logic of Computer Arithmetic, Prentice-Hall, Englewood Cliffs, N.J., 1963. 2J Greville, T.N.E., \"The Pseudoinverse of a Rectangular or Singular Matrix and its Application to the Solution of Systems of Linear Equations\", SIAM Review, Vol. 1, No. .1, January,} 1959. 3J Greville, T.N.E., \"Some Applications of the Pseudoinverse of a. Matrix\", SIAM Review, Vol. 2, No. 1, January, 1960. 4J Hanafy, A.A.R. and .Bohn, E.V., \"Rapid Digital Identifi cation of Linear Systems\", Joint Automatic Control Conference, St. Louis, Mo., August, 1971. 5cJ Ho, Y.C., \"The Method of Least Squares and Optimal Filter ing Theory\",, RAND Corp.,, Santa Monica,, Calif..,. RM-3329-PR, October, 1963. 6J Kalman, R.E., \"A New Approach to Linear Filtering and Prediction Problems\", Trans. ASME, Series D, J. Basic Eng., Vol. 82, March, 1960, pp. 35-45. 7J Kalman, R.E. and Bertram, J.E., \"Control System Analysis and Design Via the 'Second Method' of Liapunov, I. Continuous-Time Systems and II. Discrete-Time Systems\", Trans. ASME, Series D, J. Basic Eng., Vol. 82, June, 1960, pp. 371-400. 8J Kishi, F.H., \"On-Line Computer Control Techniques and Their Application to Re - entry Aerospace Vehicle Control\", Advances in Control Systems Theory and  Application, CT, Leondes, ed. , Vol. 1, Academic Press, New York, 1964. 9J Netravali, A.N. and de Figueiredo, R.J.P., \"On the Iden tification of Nonlinear Dynamical Systems\", IEEE Trans. Automat. Contr. , Vol. AC-16, No. 1, February, 66 1971, pp. 28-36. 10j Ogata, K. , State - Space Analysis of Control Systems, Prentice-Hall, Englewood Cliffs, N.J., 1968. 11j Penrose, R. , \"A Generalized Inverse for Matrices\", Proc. Cambridge Philos. Soc. , Vol. 51, 1955, pp. 406-413. 12j Penrose, R., \"On Best Approximate Solutions of Linear Matrix Equations\", Proc. Cambridge Philos. Soc. , Vol. 52, 1956. 13J Sage,-A.P., Optimum Systems Control, Prentice-Hall, Englewood Cliffs, N.J., 1968. 14J Sage, A.P. and Melsa, J.L., Estimation Theory with Ap  plications to Communications and Control, McGraw-Hill, New York, 1971. 15j Sinha, N.K. and Pille, W. , \"Online System Identification Using Matrix Pseudoinverse\", Electronics Letters, Vol. 6, No. 15, July 23, 1970, pp. 453, 454. 16j Van Trees, H.L., Detection, Estimation and Modulation  Theory, Part I, John Wiley and Sons, New York, 1968. 17j Young, P.C., \"Applying Parameter Estimation to Dynamic Systems-Part I\", Control Engineering, Vol. 16, No. 10, October, 1969, pp. 119-125. 68 APPENDIX Program-Equivalents of Symbols Appearing in the Text Programs ' Text A B C ^k \u2022 CSQU c. c, \u2014k \u2014k CTHR \u2022 threshold for c^Cj H hu * \u2014k' \u2014k I i J 3 K0 kMAX P Q ^k R r, n RS rs, n(m+n) S s, m+n SS 2 2 s , (m+n) V sk 69 Memory-Allocation for Identification Programs Labels apply to main-program assembly only. Locations avail able for user-written programs are marked by asterisks. Locations Label:Content 0000-0001 .0.0.02 0003 0004-0007 0010-0031 0040-0043 0044-0277 Use 0300 KEEP 0301 SAVE 0302 AMAT 0303 AMAT0 0304 BMAT 0305 BMAT0 0306 0307 ZERO - 0 0310 0 0311 ONE 1 0312 R 0313 S 0314 RS 0315 SS 0316 I 0317 J 0320 0321 K0 0322 L 0323 L0 0324 M 0325 N 0326 N0 0327 0330 A \u2022 500 0331 P 700 0332 Q: 1100 0333 TEMPI: 1300 0334 TEMP 2 1500 0335 TEMP 3 1700 0336 H 2100 0337 C : 2120 starting address of main program required by floating-pt. interpreter required by floating-pt. interpreter 5K &gt;pointers, indicators } floating-point zero one no. of columns in parameter array (r) number of rows in parameter array (s) product rs product ss indicator counter \u2022indicators and counters &gt;indirect matrix addresses 70 0340 B: : 2140 0341 CSQU : 2160 0342 CTHR : 040420 0343 0 0344 2700 0345 3054 0346 v' : 040420 0347 0 0350 MXADD: 2230 0351 MXSUB' 2251 035.2 ..MXMP.Y' .2.272 0353 MXDIV: 2355 0354 MXTR 2374 0355 DATRD' 2440 0356 DATPN' 2473 0357 DATRC 2533 0360 DIGIT 2572 0361 DATWR 2600 0362 WRITE 2646 0363 INIT 0364 MEAS 0365 DAT IN 0366 DTOUT 0367 0370 STR1 2170 0371 \u2022 STR2 2174 0372 STR3 2200 0373 STR4 2204 0374 STR5 . 2210 0375 \u2022 STR6 ' 2214 0376 STR7 2220 0377 STR8 2224 0400-0477 0500-2161 2162-2167 2170-2227 2230-2431 2432-2437 2440-2666 2667 2700-3403 BEGIN 3404-5577 5600-6577 ^\u2022indirect matrix addresses \"^threshold for \u00A3k\u00A3k, initially 1 contains starting adr. of main prog, contains starting adr. of calc'ns ^measurement variance, initially floaded as floating-point \"1\" &gt;indirect subroutine addresses indirect addresses for user's subroutines indirect addresses for teleprinter message strings floating-pt. interpreter work area matrix storage area (see 330-341) ( teleprinter message strings matrix arith. subr. (see 350-354) -x I/O subroutines (see 355-362) main program basic floating-pt. interpreter 71 Instructions for Using the Identification Program-Package First load the program tapes in the following order: 1. Nova Basic Floating-Point Interpreter 2. Data-supply subroutines (INIT, MEAS,DATIN, DTOUT) 3. Identification program-package The program is self-starting and will,begin by printing certain questions which are to be answered by typing numbers into the teletypewriter. Each number will be required in ei ther fixed- or floating-point format. In the case of fixed-point only one decimal digit will be accepted, while floating point format can be any string of characters in the following order: ...1. A .+ or - sign (optional-) 2. A string of decimal digits (optional) 3. A decimal point (optional) 4. A string of decimal digits 5. The letter E, if there is to be an exponent 6. A + or - sign (optional) 7. One or two decimal digits denoting exponent (optional) 8. A \"space\" A character typed in error can be deleted with a \"rubout\". Examples of allowed strings are: 500_, +50_, +5.E2_, -2.05E-04_, + . 3054E-22__, -2E03_, where _ denotes a \"space\". The questions printed and explanations of the required responses are as follows: 1. \"R = \" : A fixed-point integer from 1 to 8 equal ling r, the number of columns in the parameter array, or n, the number of system outputs. 2. \"S = \" : A fixed-point integer from 1 to 8 equal ling s, the number of rows in the parameter array, or m+m, where m is the number of system inputs and n is the number of 72 system outputs* 3. \"SAMPLES?\" : A number in floating-point format equal to the number of state-samples to be taken. 4. \"COPY? \" : A fixed-point integer corresponding to one of the following instructions regarding starting values: 0: No starting values are available for A and P_. 1: The starting values now in memory are to be used. 2: Copy the starting values from memory onto paper tape. 3: Copy the starting values from memory onto the teleprinter. 4: Read the starting values from the tape in the high-speed reader and enter them into memory (tape must be one which has been produced by response 2). 5: Accept the starting values from the teletype keyboard. (Note that following this response the program will print \"PARAMS '\" after which the elements b\"f the parameter matrix should be typed in floating-point format, row by row. A carriage-return and line-feed will occur automatically after each element has been typed and an extra line-feed will occur at the end of each row. When all rows are finished, the pro gram will print \"P-MATRIX\" and the starting values of the IP-matrix elements should then be typed row by row in floating point format.) 6: Execute the user-written subroutine whose start ing address is found in location INIT = 363. 7: Accept new initial values for CTHR and V from the teletype keyboard. (The program will respond by printing \"CTHR, V:\" after which the values of CTHR and V should be typed in floating-point format one after the other. ) 5. \"READY? \" : A fixed-point integer corresponding to: 0: Return for another pass at question 4. 1: Proceed to execute the identification program. IMPORTANT:\"The last response to question 4 must be \"0\" or \"1\" 73 Instructions for Writing I/O Subroutines INIT: This subroutine is called if a \"6\" is typed in response to the question \"COPY?\" and allows the user to supply starting values with his own subroutine. Starting address should be stored in location 363. MEAS: This subroutine is called just before the recursive identification .process .is started and can be used for such tasks as rapid pre-measuring and storing of dajta. Its start ing address should be entered into location 364. DATIN: Called each time 'a new sample of the system outputs is required. Starting address should be loaded into location 365. DTOUT: Called just after a new estimate of the parameters has been calculated and useful for outputting the parameter matrix. The starting address should be loaded into location 366. The model estimated by the.program .will be ..(.see Chapter V.) xk(n) = 0(nxn) xk_x(n) + A(nxm) uk_1(m) = A ' (nx(n+m) ) h_k (n+m) where *k = The user' s data-supply subroutines should store measured val ues of x_k and u_k in the hr-vector locations, which begin at the address found in location H= 336. xk is stored first and then u_k, each element to be written in 32-bit hexadecimal floating point format occupying 2 consecutive locations as provided by the Nova instruction FFLO. The maximum number of elements in -k-l A = A 74 h is 8, The estimated parameter array A will be left row by row starting at the address contained in location A = 330, each element in floating-point format and occupying 2 consecutive locationso Output via the teleprinter or paper-tape punch can be achieved using the subroutines which are addressed indi rectly through locations DATWR = 361 and DATPN = 356&lt;&gt; Values of the variance term s^ for a maximum-likelihood filter can be entered in floating-point format into locations V= 346 and V+1 = 347. The total number of state samples measured or estimation cycles performed is controlled by the user ' s subroutines, using location K = 320 as a counter. The initial count, which is typed in response to the question \"SAMPLES?\", is found in location K0 = 321. Dimension parameters typed in response to \"R = \" and \"S = \" and the locations where they are stored are: R = 312 r, n S = 313 s , m + n RS = 314 rs, n(m+n) SS =-315 s2, (m+n)2 Example: The following set of programs are examples of the sub routines, MEAS, DATIN and DTOUT, required to make and store a rapid set of state-samples of a continuous system via an A/D converter with multiplexed inputs, and under control of an external real-time clock. The stored data is to be proc essed one row at a time by the identification program, after which the parameter estimate is to be printed, along with a 75 warning in the case of a minimum-norm composition if an in sufficient number of linearly independent measurements were available for the parameters to be observable. ; DATA-SUPPLY SUBROUTINES FOR A/D CONVERTER R = 312 S = 313 I = ...3.1.6 J 317 K = . 320 K0 321 N 325 A = 330 H 336 BEGIN = 344 START 345 MXTR 354 DATWR = 361 WRITE = 362 STR5 374 cLOC 364 ME AS DAT IN DTOUT MAX: STORE: MEAS: SMPL: .LOC 3410 111116 123040 115105 101123 2040 3537 LDA 3, STA 3, LDA 3, STA 3, SUB 2, ADD 3, DSZ N JMP LDA R N K0 K 2 2 .-2 MAX SUBZ# 2, 3, SNC JMP @BEGIN LDA 3, STORE STA 3, 21 ; SUBO 2, 2 J LDA 3, R STA 3, N ; STRING 11: \"INS MEAS\" ; MAX NO OF STOR LOC AVAILABLE ; IND ADR FOR 1ST STOR LOC N R ; PRESET SAMPLE COUNTER ; CLEAR AC2 AC2 = KR AC3 = MAX SKIP IF KR NOT EXCEED MAX RESTART PRESET LOC POINTER AC 2 = 0 RESET MEASUREMENT COUNTER 76 CRRNT: DOAC 2, 44 NIOS 63 SKPDN 63 JMP .-1 NIOC 51 NIOS 51 SKPDN 51 JMP.o-l DIA 0, 51 STA 0, @21 NIOP 44 DSZ . N JMP CRRNT DSZ K JMP SMPL LDA 3, K0 STA 3, K ISZ K LDA 3, STORE STA 3, 21 JMP \u00A9START ; SET MUX CHANNEL TO 0 ; ENABLE HARDWARE CLOCK WAIT FOR CLOCK CLEAR A/D START A/D WAIT FOR A/D GET RESULT STORE RESULT INC MUX CHANNEL SKIP IF DONE CURRENT .SAMPLE SKIP IF DONE ALL SAMPLES PRESET SAMPLE COUNTER DATIN: STA 3, RETURN DSZ K JMP .+2 JMP 'OUT LDA 2, H LDA 3, R STA 3, N CMPNT: SUB 0, 0 LDA 1, @21 MOVL# 1, 1, SZC COM 0, 0 STA 0, 0, 2 STA 1, 1, 2 FETR FFLO 0, 2 FIC2 FEXT DSZ N JMP CMPNT JMP \u00A9RETURN RETURN: 0 STR11: 3410 OUT: LDA 3, I MOV 3, 3, SZR JMP PRINT LDA 3, J SUBZL 1, 1 ADCZ# 1, 3, SNC JMP .+3 LDA 2, STR11 JSR @WRITE ; STORE RETURN ADR ; SKIP IF NO MORE DATA ; PRESET LOC POINTER ; AC3 = R PRESET MEAS COUNTER AC0 = 0 GET DATA WORD SKIP IF NON-NEGATIVE AC0 = 177777 ; STORE DATUM IN H ; CONVERT TO FP ; INC LOC POINTER ; SKIP IF HAVE ALL COMPS : RETURN SKIP IF I = 0 AC 3 = J AC1 = 1 SKIP IF J GREATER THAN 1 ; TYPE \"INS MEAS\" 77 PRINT: LDA 2, STR5 JSR \u00A9WRITE LDA 0, S LDA 1, R LDA 2, A JSR @ DATWR JMP \u00A9BEGIN DTOUT: JMP. 0, 3 o END ; TYPE \"PARAMS PRINT PARAMETERS RESTART MAIN PROG RETURN INTERFACE FOR REAL-TIME CLOCK DS1 DS2 DS3 DS4 DS5 DS6 On STRT 3&gt; CLOCK SELD 78 Identification Programs for Nova Computer On the following pages are found the assembler listings of the basic identification programs for the Nova computer. J MATRIX ARITHMETIC SUBROUTINES ; REQUIRE BASIC FLOATING POINT INTERPRETER 79 000300 KEEP \u2022 = 300 000301 SAVE = 301 000302 AMAT 302 000303 AMAT0 = 303 000304 BMAT = 304 000305 BMAT0 = 305 000322 L 322 000323 L0 = 323 00032 4 M = 324 00032 5 N - 325 000326 N0 = 326 00030 7 \u2022 LOC 307 00 307 000000 2ER0: 0 00310 000000 0 000350 *LOC 350 00350 0022 30 MXADD 00351 002251 MXSUB 00352 00227 2 MXMPY 00353 002 35 5 MXDI V 00354 002 37 4 MXTR 002230 . LOC 2230 i ThESE-ARE THE PAGE-0 ; ADDRESSES IN WHICH THE f STARTING ADDRESSES OF i THE SUBROUTINES CAN BE ; FOUND. ; SUBR TO ADD TWO MATRICES (C = A + B) ENTER WITH: LOC N = NO. OF ELTS IN EACH MATRIX AC0 = ADR OF A AC1 = ADR OF B AC2 = ADR OF C 02230 0 40302 MXADD: STA 0, AMAT 02231 044304 STA 1* BMAT 02232 054301 STA 3* SAVE 02233 006004 XADD: FETR 02234 022302 FLDA 0* \u00A9AMAT 02235 026 30 4 FLDA 1* 6BMAT 02236 123000 FADD I* 0 02237 041000 FSTA Q* 0* 2 02240 104000 FIC2 02241 100000 FEXT 022 42 010302 I Si\u00A3 AMAT 02243 010302 ISE AMAT 02244 010304 IS2 BMAT 02245 01030 4 IS2 BMAT 02246 01432 5 DSii N 02247 000764 JMP XADD 02250 002301 JMP \u00A9SAVE ) PRESET A-MATRIX POINTER 1 PRESET B-MATRIX POINTER ; STORE RETURN ADR i ENTER FP MODE ; GET ELT OF A i GET ELT OF B ; ADD ELTS ; STORE RESULT INC i INC C-MATRIX POINTER I EXIT FP MODE ; INC A-MATRIX POINTER. I INC B-MATRIX POINTER I SKIP IF ALL ELTS ADDED ; ADD NEXT PAIR OF ELTS i RETURN I SUBR TO SUBTRACT ONE MATRIX FROM ANOTHER CC = A -; ENTER WITH: LOC N = NO. OF ELTS IN EACH MATRIX AC0 AC1 AC2 \u00BB ADR OF A = ADR OF B = ADR OF C 80 02251 040302 MXSUB: 022 52 0 4430 4 02253 054301 02254 006004 XSUB: 02255 022302 02256 026304 02257 122400 02260 041000 02261 104000 02262 100000 02263 010302 02264 010302 02265 010304 02266 010304 02267 01432 5 02270 000764 02271 002301 STA &amp;* AMAT STA 1* BMAT STA 3* SAVE FETR FLDA 0&gt; \u00A9AMAT FLDA \\\u00BB 0BMAT FSUB 1* 0 FSTA 0* 0&gt; 2 FIC2 FEXT IS 2 AMAT IS2 AMAT .IS2 BMAT IS2 BMAT . DS\u00A3 N JMP XSUB JMP \u00A9SAVE J PRESET A-MATRIX POINTER ; PRESET B-MATRIX POINTER i STOKE RETURN ADR i ENTER FP MODE ; GET ELT OF A t GET ELT OF B ;'SUBTRACT ELT OF B FROM ELT OFA I STORE RESULT IN C t INC C-MATRIX POINTER i EXIT FP MODE ; INC A-MATRIX POINTER ; INC B-MATRIX POINTER ; SKIP IF DONE ; DO ANOTHER SUBTRACTION I RETURN i SUBR TO MULTIPLY TWO MATRICES (C = AB) ; ENTER WI TH s LOC L = NO. OF COLUMNS IN A/ROWS I ; LOC M = NO. OF ROUS IN A i LOC N NO. OF COLUMNS IN B i AC0 ADR OF A t AC1 = ADR OF B 1 AC2 . \u2022 -= ADR OF C 02272 0 40 302 MXMPY t STA QJ&gt; AMAT 0227 3 0 4030 3 STA 0t AMAT0 I PRESET A-MATRIX POINTERS 02274 0 4430 4 STA 1*. BMAT 02275 0 4430 5 STA 1 * BMAT0 } PRESET B-MATRIX POINTERS 02276 054301 STA 3* SAVE t STORE RETURN ADR 02277 034322 LDA 3* L 02300 054323 STA 3* L0 i LOC L0 = NOo OF COLUMNS IN 02301 024325 LDA 1 \u00BB N 02302 0 44326 STA \\&gt; N0 i LOC N0 = NO. OF COLUMNS IN 02303 127000 ADD \\* 1 i AC1 = T'.vICE NO. COLS IN B 0230 4 f.\u2022\u2022(.'(. \u2022/.,\"(.: JMP XK PY + 2 02 305 C-34302 MRE'i : L DA 2&gt; 02.306 054303 STA 3* AMAT0 BEGIN NEXT ROW OF A 02307 03*30 5 LDA 3* BMAT0 02310 0 5430 4 STA 3&gt; BMAT GO TO FIRST COLUMN OF B 0231 1 034326 LDA 3* N0 02312 054325 STA 3* N RESET COLUMN-COUNTER 0231 3 000407 JMP XMPY 0231 4 034300 NRET: LDA 3* KEEP 0231 5 054304 STA 3* BMAT 023i6 010304 IS\u00A3 BMAT 02317 010304 IS2 BMAT j BEGIN NEXT COLUMN OF B 02320 034303 LDA 3* AMAT0 02321 054302 STA 3* AMAT i REPEAT SAME ROW OF A 02322 034323 XMPYS LDA 3* L0 02323 054322 STA 3* L i RESET PRODUCT-COUNTER 02324 034304 LDA 3&gt; BMAT 02325 054300 STA 3J&gt; KEEP i STORE COLUMN-POINTER 02 326 006004 FETR t ENTER FP MODE 02327 030307 FLDA 2* iiERO i iiERO CUMULATIVE SUM 02330 022302 LRET: FLDA 0* 0AMAT 1 GET ELT OF A 81 02331 026304 FLDA 1* @8MAT I GET ELT OF B 02332 120100 FMPY 1* 0 I MULTIPLY ELTS 02333 113000 FADD 0* 2 i ADD PROD TO CUMULATIVE SUM 02334 100000 FEXt t EXIT FP MODE 02335 010302 ISZ AMAT 02336 0 10302 IS2 AMAT i MOVE ALONG ROW OF A 02337 034304 LDA 3* BMAT 02340 137000 ADD \\\u00BB 3 02341 0 5 430 4 STA 3* BMAT \u2022 * MOVE DOWN COLUMN OF B 02342 006004 FETR t ENTER FP MODE 02343 0 14322 FDSZ L i SKIP IF ALL PRODUCTS DONE 023 44 000764 FJMP LRET i FORM NEXT PRODUCT 02345 051000 FSTA 2&gt; 0* 2 % STORE RESULT IN C 02346 104000 FIC2 i MOVE ALONG ROW OF C \"02347 100000 FEXT i E-X-I T FP -MODE 02350 01432 5 DSL' N t SKIP TF DONE ALL COLS OF B 02351 0007 43 JMP NRET 02352 014324 DS2 M t SKIP IF DONE ALL ROWS OF A 02353 000732 JMP MRET 02354 002 30 1 JMP G&gt;SAVE i RETURN 82 I SUBR TO DIVIDE A MATRIX BY A.SCALAR &lt;C = A/B) i ENTER WITH: LOC N = NO. OF ELTS IN A J AC0 ADR OF A ; AC1 ADR OF B ; AC2 = ADR OF C 02355 0 40302 MXDIVJ STA 0* AMAT PRESET A-MATRIX POINTER 02356 044304 STA 1* BMAT i PRESET B-MATRIX POINTER 02357 054301 STA 3* SAVE j STORE RETURN ADR 02360 006004 XDI VJ FETR 1 ENTER FP MODE 02361 022302 FLDA 0J \u00A7 AMAT \u2022 t GET ELT OF A 02362 026304 FLDA 1* 0BMAT I GET ELT OF B 02363 120200 \"F DT V 1* 0 i DIVIDE' ELT OF A 0236 4 041000 FSTA 0* 0J&gt; 2 I STORE RESULT' IN C 02365 10 400 0 FIC2 % INC C-MATRIX POINTER 02366 100000 FEXT i EXIT FP MODE 02367 010302 ISE AMAT 02370 010302 IS2 AMAT i INC A-MATRIX POINTER 02371 014325 DS'Z N i SKIP IF ALL ELTS DIVIDED 02372 000766 JMP XDIV 02373 002301 JMP \u00A9SAVE I RETURN J SUBR TO TRANSPOSE A SQUARE MATRIX (A = A*&gt; i ENTER WITH! AC1 NO.. OF -ROW'S OR COLUMNS OF A AC2 ADR OF A 0237 4 0 4432 4 MXTR % STA \\, M 02375 014324 DSif M ; LOC M = I LESS THAN NO. ROWS 02376 000402 JMP .*2 0237 7 001400 JMP 0* 3 02400 054301 STA 3* SAVE ; STORE RETURN ADR 02401- 127000 ADD 1\u00BB 1 ; AC I = TWICE NO. OF ROWS 02402 121400 INC 1 * 0 02403 101400 INC 0* 0 t AC0 = TWICE NO. ROWS + 2 0240 4 000403 JMP \u00BB+3 02 405 030303 TRRETS LDA 2&gt; AMAT0 02406 113000 ADD 0* 2 i MOVE DOWN DIAGONAL 02407 050303 STA 2* AMAT0 t STORE ELT POINTER 02410 050302 STA 2t AMAT i PRESET FIRST ELT POINTER 0241 1 034324 LDA 3* M 02412 054325 STA 3* N i PRESET COUNTER 02413 034302 XTR J LDA 3* AMAT 0241 4 137 000 ADD 1*3 02415 054302 STA 3* AMAT i SET FIRST ELT-POlNTER 02416 006004 FETR i ENTER FP MODE 02417 104000 FIC2 I SET SECOND ELT-POINTER 02420 026302 FLDA 1* 0AMAT 02421 031000 FLDA 2\u00BB 0* 2 I GET ELTS 02422 045000 FSTA 1* 0* 2 02423 052302 FSTA 2/ 6AMAT i SWAP ELTS 02 42 4 100000 FEXT EXIT FP MODE 02 425 01432 5 DS\u00A3 N i SKIP IF DONE ROW 02426 000765 JMP XTR 02427 01432 4 DS2 M i SKIP IF DONE MATRIX 02430 000755 JMP TRRET JMP @ SAVE ' i RETURN o END 7 777 AMAT 000 302 AMAT0 000 303 BMAT 000304 BMAT0 000305 KEEP 000300 L 000322 L0 00032 3 LRET 002330 M 000324 MRET 00230 5 MXADD 002230 MXDI V 002355 MXMPY 002272 MXSUB 0022.5.1 MXTR 002374 N 00032 5 N0 000326 NRET 002314 SAVE 000301 TRRET 002 40 5 X ADD 002233 XDI V 002 360 XMPY 002322 XSUB 0022 5 4 XTR 002413 HERO 000307 ; INPUT-OUTPUT SUBROUTINES FOR TTY AND PTP 85 $ REQUIRE BASIC FLOATING-POINT INTERPRETER 000040 00040 002560 00041 002637 -LOC RECV TYPE 40 i THESE ARE THE SUBROUTINES i 'TO BE USED BY FP INSTRUCTIONS ; FDFC AND FFDC* RESPECTIVELY 00355 00356 00357 00360 00361 00362 000300 000301 00032 4 000325 000326 000355 002 440 002 47 3 002533 002 57 2 002600 002646 KEEP SAVE M N N0 \u2022- 300 = 301 = 324 = 32 5 = 326 \u2022LOC 355 DATRD DATPN DATRC DIGIT DAT WR WRITE i THESE-ARE THE -PAG-E--0 J ADDRESSES IN WHICH THE t STARTING ADDRESSES OF I THE SUBROUTINES CAN BE i FOUND. 002440 \u2022LOC 2440 I SUBR TO STORE DATA FROM PAPER TAPE ; ENTER WITH: AC1 I AC2 NO- OF FLOATING-POINT DATA STARTING ADR OF STORAGE LOC 02440 127000 DATRD: ADD 1* 1 02 441 0 4432 5 STA 1* N 02442 050020 STA 2* 20 02 443 01 402 0 DS2 20 02444 054301 STA 3* SAVE 02445 060112 NIOS PTR 02446 063612 SKPDN PTR 02447 000777 JMP .-1 02450 060512 DIAS 0* PTR 02451 063612 SKPDN PTR 02452 000777 JMP .-1 02453 101005 MOV 0* 0* SNR 02454 000774 JMP .-4 02455 004405 JSR READ 02456 042020 STA 0* \u00A920 02457 014325 DSZ N 02460 000775 ' JMP .-3 02461 002301 5 JMP \u00A9SAVE 02462 064512 READ: DIAS 1* PTR 02463 063612 SKPDN PTR 02464 000777 JMP .-1 02465 125300 \u2022 MOVS 1* 1 \u2022\u2022 02466 060512 DIAS 0* PTR 02467 063612 SKPDN PTR 02470 000777 JMP .-1 02471 123000 ADD 1* 0 02472 001400 ' JMP 0* 3 ; DOUBLE AC1 t N = NO o OF DATA WORDS I SET LOC POINTER ; STORE RETURN ADR i READ A LINE FROM TAPE ; GET RESULT* READ AGAIN i SKIP IF RESULT N0N-2ER0 t GET DATA WORD J STORE DATA WORD ; SKIP IF ALL WORDS READ ; RETURN I GET RESULT* READ AGAIN i LEFT-JUSTIFY IN AC1 i GET RESULT* READ AGAIN I COMBINE HALVES i RETURN ; SUBR TO PUNCH DATA ON PAPER TAPE 86 1 ENTER WITH: AC1 i AC2 = NO. OF FLOATING-POINT DATA = STARTING ADR OF STORAGE LOC 0247 3 127000 DATPN: ADD 1* 1 ; DOUBLE AC1 02474 044325 STA 1* N ; N = NO. OF DATA WORDS 02475 050020 STA 2\u00BB 20 02476 014020 DS2 20 t SET LOC POINTER 02 47 7 '0 54301 STA 3* SAVE i STORE RETURN ADR 02500 102 400 SUB 0* 0 i 2ER0 AC0 02501 152 420 SUB2 2&gt; 2 i 2ER0 AC2* SET CARRY 02502 004421 JSR PUNCH t PUNCH A 2ER0 02 503 151103 MOVL 2* 2* SNC i COUNT OF 17 -02-504 \u202200 0-7-7-6 JMP -0-2 02505 100000 COM 0* 0 t SET AC0 02506 061113 DOAS 0* PTP i PUNCH A 377 02507 063613 SKPDN PTP 02510 00077 7 JMP \u00BB-l 0251 1 022020 LDA 0* @20 i GET DATA WORD 02512 004411 JSR PUNCH i PUNCH DATA WORD 0251 3 01432 5 DS2' N i SKIP.IF ALL WORDS PUNCHED 0251 4 00077 5 JMP .-3 02515 1 02400 SUB 0* 0 i 2ERQ AC0 0251 6 1 52 420 SUBE 2* 2 i 2ER0 AC2* SET CARRY 02517 00440 4 JSR PUNCH i PUNCH A 2ER0 02520 151103 MOVL 2* 2* SNC \u2022 * COUNT OF 17 02521 00077 6 JMP .-2 02522 002301 JMP QSAVE i RETURN 02523 105300 PUNCH % MOVS 0* 1 02524 065113 DOAS 1* PTP 02525 063613 SKPDN PTP 02526 000777 JMP .-1 02527 061113 DOAS 0* PTP 02530 063613 SKPDN PTP 02531 000777 JMP -.-1 02532 001400 JMP 0* 3 \u2022 I RIGHT-JUSTIFY FIRST HALF i PUNCH FIRST HALF i PUNCH SECOND HALF t RETURN I SUBR TO STORE DATA I ENTER WITH: AC0 J AC! i AC2 02533 04032 4 DATRC % STA 0* M 02534 044326 STA 1 * N0 02535 0 5 4301 STA 3* SAVE 02536 034326 NXTRW'8 LDA 3* N0 02537 054325 STA 3* N 02540 020503 LDA 0* LF 02541 004476 JSR TYPE 02542 020502 NXTEL? LDA 0* CR 02 543 004474 JSR TYPE 02544 020477 LDA 0* LF 02545 004472 JSR TYPE 02 546 006004 FETR 02547 124000 FDFC : 1 FROM KEYBOARD -\u2022 NO. OF ROWS OF FP DATA = NO. OF COLUMNS = STARTING ADR OF STORAGE LOC J M = NO. OF ROWS i N0 = NO. OF COLUMNS ; STORE RETURN ADR J N = NO. OF COLUMNS ; LINE-FEED i CARRIAGE-RETURN i LINE-FEED . J ENTER FLOATING-POINT MODE i ACCEPT DEC NO.* CONVERT 02550 02551 02552 02553 02554 02555 02556 02557 045000 10 4000 100000 014325 0007 6 6 014324 000760 002301 FSTA 1* 0* FIC2 FEXT DS2 N JMP NXTEL DSZ M JMP NXTRW JMP \u00A9SAVE 87 I STORE HEXADECIMAL NO. I INC STORAGE-LOC POINTER ; EXIT FP MODE J! SKIP IF HAVE ALL ELTS OF ROW I SKIP IF HAVE ALL ROWS i RETURN 02560 054300 RECVJ ' STA 3* KEEP 02561 060110 NIOS TTI 02562 063610 SKPDN TTI 02563 000777 JMP .-1 02564 060410 DIA 0* TTI 02565 024404 LDA b MASK 02566 123400. AND \\* 0 02567 004450 JSR TYPE 02570 002300 JMP \u00A9KEEP 02571 000177 MASK I 177 t STORE RETURN ADR I ENABLE KEYBOARD I WAIT FOR CHARACTER i GET CHARACTER ; AC1 = 177 i MASK TO 7 BITS i ECHO CHARACTER t RETURN 02572 054301 02573 004765 02 57 4 02 4403 -02-57-5 123 400 02576 002 301 02577 000017 I SUBR TO ACCEPT A DIGIT FROM KEYBOARD i BINARY VALUE OF DIGIT IS LEFT IN AC0 DIGITS STA 3* SAVE , JSR RECV LDA 1* DTMSK -AND \\* 0 JMP \u00A9SAVE DTMSK s 17 i STORE RETURN ADR I RETURN DIGIT IN AC0 I AC1 = 17 \u2022I -MASK TO 4 .BITS ; RETURN ; SUBR TO TYPE DATA ON TELEPRINTER 88 i ENTER WITH: AC0 = NO. OF ROWS OF FP DATA i AC1 NO. OF COLUMNS i AC2 = FIRST ADR WHERE DATA STOf 02600 0 4032 4 DATWR! STA 0* M t LOC M = NO. OF ROWS 02601 0 4432 6 STA 1J&gt; N0 i N0 = NO. OF COLUMNS 02602 054301 STA 3* SAVE i STORE RETURN ADR 02603 125112 ROW: MOVLtf Is l\u00BB S2C I. SKIP IF 2 LINES TYPED 02604 000403 JMP o-f-3 i NO LINE-FEED 02605 020436 LDA 0* LF 02606 004431 JSR TYPE t LINE-FEED 02607 02 4427 LDA 1* COLS I AC1 = -4 026T0 0-3432*6 LDA 3* lN0 0261 1 054325 STA 3* N I N = NO. OF COLUMNS 02612 020 432 LINE; LDA 0* CR 02613 004 42 4 JSR TYPE t CARRIAGE-RETURN 0261 4 020 42 7 LDA 0* LF 0261 5 00 4422 JSR TYPE i LINE-FEED 02616 020427 ELTs LDA 0* SP 02617 00 4420 JSR TYPE i SPACE 02620 00 4417 JSR TYPE SPACE 02621 00600 4 FETR i ENTER FP MODE 02622 021000 FLDA . &amp;\u00BB ids 2 i LOAD HEX FP NO. IN FAC0 02623 140000 FFDC 0 i TYPE NO. IN DEC FP 02624 10 4000 FIC2 \u2022 &gt; INC STORAGE-LOC POINTER .02625 \u2022,1.0 0.0.0.0 .FEXT * EXIT FP MODE 02626 014325 DS\u00A3 N SKIP IF DONE ALL ELTS OF 02627 000 40 4 JMP \u00AB + 4 02630 01 432 4 DS?\u00A3 M i SKIP IF DONE ALL ROWS 026 31 0007 52 JMP ROW 02632 002301 JMP \u00A7SAVE i RETURN 02633 12540 4 INC \\* 1* S2R \u2022 i SKIP IF FINISHED LINE 02634 000762 JMP ELT 02635 000755 JMP LINE 02636 177774 COLS 5 -4 02637 061 111 TYPE; DOAS 0* TTO i TYPE CHARACTER 02640 063611 SKPDN TTO 02641 000777 JMP .-1 026 42 001400 JMP 0* 3 i RETURN 02643 000012 LF: 12 02644 000015 CR: 15 02645 000040 SP: 40 02646 0 54301 02647 020775 02650 004767 02651 020772 02652 004765 02653 004764 02654 024762 02655 021000 i SUBR TO TYPE A STRING OF 8 CHARACTERS i ENTER WITH AC2 = STARTING ADR OF STRING WRITE: CHAR % STA 3* SAVE LDA Q&gt; CR JSR TYPE LDA 0* LF JSR TYPE JSR TYPE LDA \\\u00BB COLS LDA 0* 0* 2 i STORE RETURN ADR i CARRIAGE RETURN ; DOUBLE-SPACE i AC1 = -4 i GET WORD 026 56 101300 MOVS 0* 0 02657 1012O0 MOVK 0* 0 02660 004757 JSR TYPE 02661 021000 LDA 0* 0* 2 02662 0047 55 JSR TYPE 02663 151400 INC 2* 2 02664 125404 INC 1* 1* 02665 000770 JMP CHAR 02666 002301 JMP \u00A9SAVE 007777 \u00BB END 7777 i SWAP HALVES I SHIFT RIGHT I PRINT 1ST CHARACTER t GET WORD AGAIN I PRINT 2ND CHARACTER I INC LOC POINTER ; SKIP IF DONE i RETURN CHAR 002655 COLS 002636 CR 002644 DATPN 002 47 3 DATRC 002533 DATRD 002 440 DATWR 002600 DIGIT 002572 DTMSK 002577 ELT 002616 KEEP 000300 LF 002643 LINE 002612 M 00032 4 MASK 002571 N 000 32 5 N0 0'00'32o NXTEL 002542 NX TRW 002536 PUNCH 002 52 3 READ 002 462 RECV 002560 ROW 002603 SAVE 000301 SP 002645 TYPE 002637 WRITE 002646 330 1 .MAIM, ; RECURSIVE LEAST-SQUARES IDENTIFICATION ; REQUIRES! ^ASIC EP INTERPRETER ; MATRIX ARITHMETIC SUBROUTINES ; I/O SUBROUTINES FOR TTY, TAPE ; DATA,-SUPPLY PROGRAMMES ; THIS PROGRAMME USES THE MAXIMUM ' ; HO. OF SYMBOLS ALLOWED 3Y 'THE ; RELOCATA3LE ASSEMBLER. DO NOT ADD ANY MORE. 3333 13 0003 1 4 0303 15 30 03 16 0iV33 17 0 rj r: 3 21 300322 030324 333325 .'\u2022 o. \"\u2022. 7, r -&lt; u r.-1.- O ti/ 33035 1 033352 300333 303354 300355 i', o .&gt; \u2022&gt; r. n \u00BB v n O i \u2022 i r i O J I 333364 333365 ^, r\\ x r c \u2022q S RS I J K3 L M M MX ADR MXSU3 y. y V] p V MXDIV MXTR DATRD. DATPM DATRC DI GIT DATWR I M IT ME AS DAT I!-! DTOUT 3 12 3 13 3 1/1 3 15 3 1 6 3 17 32 1 322 324 325 3 52 35 I 3 52 353 3 54 355 3 5 6 357 3 6 3 351 3 62 3G3 3 64 3 65 366 00002 002344 .LOC 2 JMP 3344 033037 M n \u2022 *y f ' - ^ /. r* '.*. K&gt; \u00BB; v.' ( L&gt; V..1 \u00A3. 3303 1 1 3031 1 2C0361 ONE ? .LOC 7 403 .LOC 3 1 1 1 ;. WORK AREA FOR F? KIT 0SS'333 .LOC 33 3 0333 3 ;7 &lt;T; *\u2022&gt; C O -V \u00ABfi 7 &lt;J At 533 0333 1 &lt;-\\ r.\\ 7 ^ '-. t.- &lt;L- o / i ' o . 700 30332 331 10:\" O \u2022 1 133 33333 32 133C: TEMP 1: 1333 3333 4 3 ?! 5 3 3 TEMP2: 1533 30335 .3-0 172 0 TEMP3: 17 r 33336 K: 2133 33337 3 02120 C: 2 1 23 033 4 0 3 32\\&amp;? P \u2022 2 14 0. 0334 l 3321 S3 CSRU- 216.\" i. i.' r..' 0^:3. .LOG 3 42 ; MATRIX ADDRESSES 0032 00342 00343 003 4 4 0034? 00346 00347 33370 0037 1 00372 0037 3 0037 4 00375 00376 00377 MAIM 04 0 4 20 000000 00034 4 002700 00305 4 0404^0 000000 000370 00 2170 002174 0 02200 002204 f.;'' O 9 J rf 002214 0 0 2224 CTHR: STR 1 SIR 2 SIR 3 SIR 4 SIR 5 SIR 6 STR 7 0 40420 ry. .LOC 34 4 BEGIN' START 04 04 20 0 .LOC 370 2170 2174 2 20 0 2204 2210 22 14 2220 999 h t ADDRESSES OF STRINGS 002170 .LOC .2 170 \u2022 02170 122040 1220 4 0 \u2022 0217 1 075 042 \u2022 075040 02172 f&gt; /, r? f\u00BB /; P 02173 042040 0 4 22 4;; 02174 1 2 3 0 4 2 123040 \u2022 02173 275::A0 075043 02176 0 4 2 0 4 2 if, /, n. rt /j 02177 2 4 2 2 4 2 043042 32223 123 12 1 12310 1 \u2022 0220 1 115122 115120 0220 2 1 1-4 125 .1 1 4 1 0.5 02203 12307 7 123277 r&gt; o o/ 123 1 17 103 1 17 \u2022 0??05 !20!3 1 12013 1 02205 27 7 24,;- 077042 02207 0 4 3 2 4 0 0 4 00/10 0221 0 120 121 120 10 1 \u2022 0221 1 122121 12210 1 0?? I o 115 123 1 15123 02213 0 422 4 0 04 0O40 02214 120255 120255 \u2022 02213 115 12 1 11510 1 02216 124 122 124 122 02217 111132 111130 ['9 99 (/, 122135 122105 * 1 02221 10 1104 101 104 02222 131077 131377 02223 04^040 043040 0 2224 103124 103 124 \u00BB 02225 110 122 110 122 02226 0 5 4 2 4 0543 4 2 02227 126 27.2 126072 002702 .LOC 2723 \u00AB 1 02700 0.06005 BEGIN: FIMI 0270 1 032370 LDA 2, STR1 02702 2263 62 JSR 3VRITE \u2022 J 02703 006360 JSR 'DDI GIT \u2022 02704 0 4 0312 STA 0, R MESSAGE STRINGS IM ASC : STRING I: \"R = ; STRING 2: \"S ; STRING 3: \"SAMPLES? t STRING 4: \"CO.^Y? ; STRING 5: \"PARA^S \" \u2022 STRING 6: \"R-MATRI: ; STRING .7: \"READY? { STRING 3: \"CTHR, V: ; PROGRAMME \"-EGIMS 0003 , (.270 3 02706 02707 027 10 027 1 1 027 12 027 13 027 14 027 15 G27 IS 027 17 02723 02721 22722 02723 0 272/! 02725 0 2726 02727 227 3 0 0273 1 02732 02733 0 2734 0273 5 02735 02737 K274K 027 4 1 0 27/! 2 02743 0274/1 02745 02745 027/7 0 27 5 0 027 5 1 02752 K2753 027 5 ^ 0 2755 02756 0.275 7 02760 0276 I 02762 027 63 0276/-02765 02766 0 2767 02770 0277 1 0277 2 B2773 027 7/J 0 2775 02776 02777 MAIN 03037 1 006362 006360 0 4 03 13 0/0325 126/03 107000 014325 n r.\\ in c C- LJ I / o 0 A A 3 1 5 03 43 12 0 5 43 25 1 26 4 00 107000 0 14325 0 0 0 7 7 6 0 4 4 3 1 A 030372 '\u2022\u2022 c y c o 102520 126520 030333 0 0 600 4 075000 100 00 0 035 0 0 1 93 054 1 0303 7 3 006362 006360 0403 16 10 100 5 0 004 7 6 0403 17 014317 r&gt; r\\ &lt; \\ i. n o 1 . \u00AB'.\u2022 i. - &amp; 000472 0 14.3 17 000407 024.3 14 030330 006356 0243 15 03033 1 006355 0143 17 0004. 15 030374 006352 0203 13 02/-3 12 030330 00 636 1 030375 006362 0203 13 0243 13 COPY: OPT op jo. 0PT3; JSP JSP STA ST A SUP ADO DSZ JMP STA LDA STA SUB ADD DSZ \u00ABj 11 &lt; S T ^ LDA JSP SUP7. SUPZ LDA JSR LDA FETP FFIX FEXT LDA STA LDA JSP JSP STA MOV JMP STA DSZ JMP JMP 0S7 JMP LDA LDA JSR 2 , S T P 2 f:&gt;!&gt;.! ? I T~ 001 GIT 3, M 2 1, 3, 3, 1, SS P M 1 1 _ 2 1 , PS 2, STP3-SVRITE L 1 , 1 TEMP 1 o '-\u00BB f!&gt;^ T 3 C o TEMP 1 3, K0 2, SIP 4 PWRITE f-D I 01 T I ..., iJ, PEADY a I \u00BB , J J OPT 2 PEADY J 0PT3 PS A \u00A9DATPN SOP 1 9 2, LDA 1, SS LDA JSP DSZ I WD t_M ! i LDA JSP LDA LDA LDA JSP LDA JSP LDA LDA 2, P 50ATPM J OPT A 2, STP5 (SVP ITE 0, S 1, R 2, A ; CDATVR 2,' STP6 \u00A9WRITE n. c c. , t&gt; 1 , s ; PRINT \"5 \u2022= ? ; GET S j SS = SOU APE OF S ; LOC PS - PRODUCT OF P AND 3 ; PRINT \"SAMPLES?\" ; GET KO \u2022 KO = NO. OF RAPID SAMPLES ; PR IMT \"COPY?\" : GET I ; SKIP IF I - 0 ; MO ST ART I MO VALUES \u2022 SAME STARTING VALUES ; PAPER-TAPE COPY ; TELETYPE COPY 0004 .MAIM 03033 03033 1 0300 1 00636 1 03002 01/i317 OPT/;: 03003 000^07 03004 0243 14 03005 030330 03005 006355 .03007 024315 03010 030331 030 11 00S355 03012 0 M317 0?T5? 030 13 0004 15 03014 030374 03015 005362 030 1-6 -0203 13 03017 024312 .03020 030330 030 21 006357 03022 03037 5 03023 006362 03024 0203 13 03025 024313 03026 030331 03027 006357 03030 014317 0?T6: 0303 1 000402 03032 006363 . 03033 0 14317 0PT7t 03034 0004 | ! . 03035 030377 0 3 0' 3 6 0 6 3 6 2 03037 006004 03040 120000 03041 124000 03042 040342 03043 044346 03044 100000 0304.5 030376 HEADY! 03046 006362 03047 006360 03050 10 1005 0305 1 00067 1 03052 002364 03053 003260 - 94 LDA 2, P JSP \u2022ODMWP. . DSZ. J e STARTING VALUES FROM JMP OPT5 LDA 1, PS LDA 2, A JSP \u00A9DATRD LDA 1, SS LDA 2, P . * JSP \u00A9DATRD-DSZ J \u2022 STARTING VALUES FROM JMP 0PT6 LDA 2, STR5 JSR \u00A9WRITE \u2022 PR I NT \"P ARAMS\" IDA 0.., .S LDA 1, R LDA 2, A JSR \u00A9DATRC LDA 2, STR6 JSR \u00A9WRITE \u2022 PRINT \"P-MATRIX\" LDA 0, S LDA 1, S LDA. 2, ? JSR \u00A9DATRC DSZ J \u2022 STARTING VALUES SU??L JMP .+2 . \u2022 S 3Y USER SURROUTI v \u2022? JSR \u00A910 IT DSZ J \u2022 CHANGE CTHR AMD V JMP READY LDA 2, STR--JSR \u00A9WRITE PRINT- \"CTHR , V: \" EETR FDFC 0 \u2022 GET CTHR FDFC 1 \u2022 ; GET V FSTA 0, CTHR \u2022 STORE CTHR FSTA 1, V \u2022 STORE V FEXT LDA 2, STR7 JSR \u00A9WRITE \u2022 PRINT \"READY?\" JSR \u00A9DIGIT MOV 3, 0, SMR JMP COPY JMP \u00A9MEAS \u2022 * USER PROGRAMME SEQ2 0005 030 5 4 0 3055 03056 03057 03060 030 6 1 03062 03063 030 6 4 03065 03066 03067 03070 0307 1 03072 03073 0 307 4 0307 5 0307 6 03077 0 3 1 0 0 . 0310 1 03 104 03 105 03 106 03107 03 1 10 \u00A33 1 1 1 0 3112 03 1 13 03114. 031 15 031 15 031 17 03 120 03 12 1 03 122 03123 03124 03 125 03 126 MAIM 95 006365 START \u2022 JSR \u00A9DAT IN . \u00AB 0 3 43 16 LDA 3, I 175004 MOV 3, 3, SZR \u00AB 0 02774 JMP \u00A9START-1 034313 LDA 3, S 0543 17 STA 3, J \u2022 10 2400 SUP O 9 0 30333 LDA 2, TEMPI 034330 L D A 3, A 0 4 14^0 STA 0 , 0, 3 1754.00 IMC 3, 3 1724 14 SU3\u00A3 3, 2, SZR O0O773 JMP .-3 020332 LOO?: LDA 0 0 \u00AB 024336 LDA 1 I H 050337 LDA ?., C 0343 13 LDA 3, S 054322 STA 3, L 054324 STA 3, M 0343 1 1 L D \\ 3, ONE 05 4 3 25 STA 3, N 006332 JSR \u00A9MXMPY 020335 L D A 0, H \u2022 0 24337 LDA 1 , C 030337 L D A 2, C 0343 13 LDA 3, 5 034325 STA 3, N 00633 1 JSR \u00A9'&lt;*SU3 ;203 37 LDA 0, C \u2022 024337 LDA 1,0 0 3 034 1 LDA 2, CSO.U 0 3 43 13 LD A 3, S 05 4 522 STA 3, L 034311 l?\\ 3, ONE 0 5 4 3 24 . STA 3, M 0 54325 STA 3, N 006352 JSR ,j;yjv,v. p y O 0 C n r? /. rr TO 020342 FLDA 0, CTHR 026341 FLDA, 1 , \u00A9C32U ^0403 - D ; 1. OH ; 2. C .r H - QH ; 5. C \\C C TOO SMALL? rSUB \u2022JMP TXT .+3 1, FSLE 03127 000725 JMP s TART 03 150 100000 FEXT 0315 1 0:.' T, % 7 c r n 1 . LDA n. , c. 05152 024541 LDA 1 , C32U 03 133 050540 LDA c \u00BB 3 03 13 4 0343 15 LDA' 5 . s 03133 0 5 45 25 STA 3 \u00BB l\" 03 156 006555 JSR \u00A9 \u2022IX DIV 03157 020555 LDA u 05140 0243 40 LDA 1 T&gt; 05141 0 30555 LDA T7MD , J U \u2022. 1 1 05142 034311 LDA 5 , ONE-05 145 054322 STA 3 \u00BB L 05 144 0345 15 LDA 5 \u00BB 5 05145 054524 STA 3 \u2022 t i 05 146 054325 3 : 5. 3 r C/C *C ; 6. H3V 3036 . :\u2022\u2022 A \\ M 03M7 006352 05 153 02035 1 03 15 1 024.335 03 !52 050534 03 153 3545 13 03 154 054322 03 155 054324 03 155 054325 03 1 57 0 063 5 2 03 160 02035 1 0316 1 0:245 3 4 03 162 03033 1 03 163 0 3 4 3 15 03 164 054325 03 165 00635 1 03 166 3243 13 03 167 030354 03 170 006354 03 17 1 020354 03 172 324335 03 175 r&gt; 7 7 7 ^ r.- O .,\u2022 O O J 03 174 0543 13 03175 354322 03176 334324 03 177 354335 03220 336352 0320 1 323340 03202 324342 03205 0.3.233-5 03204 3343 1 1 05205 254522 03206 0345 15 03207 0 5 43 24 ,03210 054525 0321 1 217\u2022\u2022 6 3 5 9 03212 0.2053 1 03213 024354 03214 03033 1 03215 0343 1 5 03216 05 4325 0 3217 006351 03220 020351 0 3221 024355 03222 35055 1 03225 0545 15 03224 05 4525 05223 &lt;-\u2022 .'\u2022 n t, i:. O 0 -&gt; *' 03226 223331 03227 024355 03230 333551 0325 1 0543 15 03232 054325 05255 006353 03234 020340 05255 024357 03236 333333 03237 0545 1 1 03240 3 5 4 7, 2 2 0324 1 0343 13 JSR 3XXMPY LDA D 1 LDA 1, T r v P 1 LDA 0 t TEXP2 LDA 3, S STA 3, L STA 3, M STA 3, M :\\ JSR :XPY LDA &gt; P LDA 1, TEXP2 LDA 2, p LDA 3, SS STA 3, V JSR 33 y. \u2022.SL'R LDA \u2022s LDA 0 - * TF.^1?2 JSR . T -:&gt; . ~ &gt; \\ LDA f \" \u00BB T E Y P 2 LDA 1, TEX 0 1 LDA 0 :- t T p v) 0 7. LDA .3, S STA 3, L STA 3, M STA 3, M ' JSR m y-:MPY LDA '\u20AC\u00BB B LDA -&gt; LDA 0 TEXP 1. LDA 3, 0 M E STA 3, L LDA 3, S STA 3, X STA 3, ,M JSR ry*\\ ;X pv LDA 0,' ? LDA 1, TEXP2 LDA 2, P LDA .3, SS STA 3, N JSR 3XX.SUB LDA \"\u00BB P LDA TEXP 3 LDA. 2, P LDA 3, SS \u2022 STA 3, ri JSR PXXADD LDA u, ? LDA 1, TEXPl LDA 2, 0 LDA 3, SS STA 3, M JSR S'XX ADD LDA B LDA 1, C LDA 0 TEXPl LDA 3, ONE STA 7 \u2022 \u00BB L LDA 3, S 96 j .7, PH3' ; 8. ? PHB' ; 3. (3H'P).(H3') ; 10. BB' \u2022 11. (P - PHB ' ) - RH 'P \u2022 12. (P-PHB'-3H'P)+BH'PHB ' !3. ?-(P-PHB '-RH'*P +BH *PHB ' )+B3 ' \u2022 1/ pr \u00BB 007 . MAIN 054524 STA 5, M 03243 054325 STA 3, M 03 244 006552 JSR \u00A9MXMPY \u00A33245 020532 LDA 0 0 03246 024533 LDA 1 , TEMPT 0 324 7 030352 LDA 2, Q 0 3250 034315 LDA 3, SS 0325 1 054325 STA 3, M 03252 006350 JSR \u00A9MX ADD 03253 004 465 JSR EST 03254 0 14317 DSZ J 0325 5 0 005 14 JMP LOO? 03256 0343 1 1 LDA 3, 0 E 03257 0 5 4316 STA 3, I 03260 02033 1 SEQ2: LDA 0, P 0326 I 0 24 35-6 \u2022 LDA 1, H 03262 0 30555 LDA 2 TEM ? 1 0326 3 0 5 45 15 LDA i\\ s 03264 054322 STA 3, L 03265 0 5 43 24 STA 3, M 03266 0543 11 LDA 3, ONE , 03267 054325 STA 3, N 03270 006352 JSR \u00A9MX^ 0327 1 r- O n 7, 7 c LDA 0, H 03 272 024333 1. D A 1, TEMPI 0 3273 030334 LDA 2, TEMP2 03274 0343 13 LDA 3, S 0 3275 0 5 4322 STA 3, L 05276 0345 1 i .LDA 3, 0 N 7. 03277 054324 STA 3, M 03500 05 4325 S T A 3, M 0330 ! 006352 JSR 05502 006004 FETc: 03505 0203 46 FLDA 0, V 05 50 44 025354 FLDA 1 , \u00A9TEMP2 05505 107000 . F ADD 0, 1 05306 04 555 4 FSTA 1, \u00A9TEMP2 05507 i r&gt; r* n. r? 1 * \u2022 ii- \u00A3\u2022 L. XJ FEXT 035 10 0 20555 LDA 0, TEMP 1 053 1 1 024354 LDA 1, TEMP2 05512 050540 LDA 2, B 035 15 03 43 13 LDA 3, S 055 14 \u00A354525 STA 3, M 055 15 006355 JSR \u00A9MXDIV 05516 0 20533 LDA 8, TEMPI 055 17 0245 4 0 LDA 1, B 05520 030355 LDA 2, TEMP5 05521 0545 1 I LDA 5, ONE 05522 05 4522 STA 5, L 05323 05 45 15 LDA 5, S 05324 0 5 45 24 STA 5, M 03525 054525 STA * M ^ &gt; \" 05326 006552 JSR pMXMpY 35527 0 2055 1 LDA 0, P 0 3350 024555 LDA 1, TEMP5 0533 1 05053 1 LDA 2, P 05332 0345 15 LD A 3, SS 03533 05 4325 STA 3, '\u2022! 0355 4 0 0 655 1 JSR \u00A9MXSUB 97 15. 0 r Q +, 16. I I1 17. H *PH 13. V + H ' PH 19.3 = PH/CV+M 'PH ) 20. (PH)B' 2 1. P = P - PH3' .MAI N 98 03335 00 4 4 03 JSR EST 03336 000722 JMP SEQ2 03337 r&gt; r&gt; (* r\\ n (f v.- Oiv o o z-' 0 03340 054777 EST; STA 5 , .-1 0334 1 LDA 0, H 0334 2 .0 24 330 LDA 1 \u00BB A 0334 3 030334 LDA 2 TEMP2 0334 4 03 43 13 LDA 5\u00BB S 03345 05 4 322 STA 5S L 33346 03 43 1 1 . LDA 5 , ONE 033 4 7 05 4324 STA 5, M 03350 0343 12 LDA 3, 0335 1 05 4 3 25 STA 3 t N 03352 006352 JSP XMPY 03353 006365 ' JSR \u00A9D ATI M 03354 020336 LDA - r?i H 03355 02433 4 LDA 1 ' TEMP 2 03356 030335 LDA o TEMP 3 03357 0343 12 LDA 3,' R 03360 0 5 4325 STA M .1 0336 1 006351 JSR \u00A9MXSUP 03362 02034.;; LDA ?\u00BB 03363 024 33 5 LDA 1 \\ TEMP 3 0336/J 030333 LDA 2, TEMPI 0 3365 0343 11 . . LDA 3, 0ME 03366 0 5 43 22 . STA 3, L 03367 0 3 4 3 13 LDA 3, S 03370 054324 STA 3, M 0337 1 0-345 12 LDA 3 . . p 03372 054325 STA 3 , N 03373 0 0 63 52 JSR \u00A9M XMPY 0337 4 020330 LDA 0 , 033 7 5 0 24355 LDA 1 , TEMP 1 03375 03 0350 LDA 2, A 03377 05 45 14 LDA 3, RS 03 4 00 05 4525 STA 3 , 111 0340 1 006550 JSR \u00A9MX ADD 034 02 r \"c 7 c c JSR \u00A9DT0UT 034 03 002754 JMP P~ ST- 1 ; 22. H VA ; 23. Z' - H ' A ; 24. 3 (Z '-H ' A) ; 25. A= A+3 (Z ' -1! ' A) ; USER PR0GR AM MS VJ YJ C / v z.; .END 27 0 0 000D '.MAIN A 000330 B 0003^0 BEG I M 002700 C 000337 COPY 0 0 274 2 CSO'J 000341 CTHR 000342 DAT IN 000305 DATPN r* n r\\ x \u00A3\\ c . t, &lt;.:- &lt;. o j o DATRC 000357 DATRD 000355 DATWR &lt;\u2022 i. &lt;.\u2022 \u2022_) 1 DIGIT / r&gt; r' i r r* i. -v &gt;\u2022\u2022 o O t; DTOUT 000366 EST \u2022 003342 H 000336 I 0003 ! 6 INIT 000363 J 0003 17 K0 L 000322 LOOP 00307 1 M 000324 MEAS 0C0 354 MX ADD 000350 MX DIV 030353 MXMPv 000352 MXSUB 000351 XX TH 00 0-3-51 N 000325 ONE . 0003 1 1 OPT1 002751 0RT2 002754 OPT 3 002764 0PT4 003002 OPT 5 0030 12 OPT6 003030 OPT 7 003033 o I 00033 1 0 000332 R 00 0312 READY 003045 RS 0 003 U S 000313 SEQ1 003 13 1 SEQ2 003260 SS 0003 15 START 003054 STR 1 000370 STR2 00037 1 STR 3 000372 STR 4 0003 73 STR5 000374 STRS 000375 STR 7 000376 STR8 000377 TEMPI 000333 TEMP 2 000334 TEMP 3 000335 06 10 .MAIN 100 V 0003^6-WRITE 00036R "@en
<https://open.library.ubc.ca/cIRcle/collections/ubctheses/831/items/1.0101755>