|
|
// EXEC FORT1CLG,REGION=512K 00000010 //* 00000020 /*JOBPARM COPIES=4 00000030 C 00000040 C PROGRAM GFLS: 00000050 C EQUATION NUMBERS IN COMMENT STATEMENTS REFER TO 00000060 C ''FLEXIBLE LEAST SQUARES FOR APPROXIMATELY LINEAR SYSTEMS'' 00000070 C BY R. KALABA AND L. TESFATSION, IEEE TRANSACTIONS ON SYSTEMS, 00000080 C MAN, AND CYBERNETICS 20 (1990), PAGES 978-989. 00000090 C LAST UPDATED: 14 JUNE 1994 00000100 C 00000110 IMPLICIT REAL*8(A-H,O-Z) 00000120 INTEGER T,TCAP,TCAP1 00000130 REAL*8 M 00000140 C 00000150 C THIS PROGRAM IS CURRENTLY DIMENSIONED FOR A MAXIMUM OF TCAP=110 00000160 C OBSERVATION VECTORS Y OF MAXIMUM DIMENSION MOBS = 15 WITH STATE 00000170 C VECTORS X OF MAXIMUM DIMENSION N = 15. 00000180 C 00000190 DIMENSION QO(15,15),PO(15,15),RO(15,15),QZERO(15,15) 00000200 DIMENSION PZERO(15,15),RZERO(15,15) 00000210 DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15) 00000220 DIMENSION M(15,15),Y(15,15),TRUEX(15,110),YY(15,110) 00000230 DIMENSION HT(15,15),U(15,15),C(15,15),W(15,15),V(15,15) 00000240 DIMENSION E(15,15),Z(15,15),G(15,15),QNEW(15,15),PNEW(15,15) 00000250 DIMENSION S(15,15),RNEW(15,15),XTCAP(15,15),X(15,110) 00000260 DIMENSION AA(15,15),BB(15,15),CC(15,15),DD(15,15),EE(15,15) 00000270 DIMENSION FF(15,15),HH(15,15),OO(15,15),PP(15,15),QQ(15,15) 00000280 DIMENSION RR(15,15),TT(15,15) 00000290 C 00000300 C ADDITIONAL ARRAYS IF SMOOTHED ESTIMATES ARE TO BE CALCULATED 00000310 C FOR INTERMEDIATE X VALUES (I.E., IF IFLAGS IS SET AT 1) 00000320 C 00000330 DIMENSION GG(15,15,110),SS(15,110) 00000340 C 00000350 C THE FOLLOWING SUBROUTINE INITIALIZES THE PENALTY WEIGHT AMU, 00000360 C THE NUMBER OF OBSERVATIONS TCAP, THE STATE VECTOR DIMENSION 00000370 C N, THE OBSERVATION VECTOR DIMENSION MOBS, AND THE INITIAL COST 00000380 C FUNCTION CHARACTERISTICS QZERO, PZERO, AND RZERO. IT ALSO SETS 00000390 C THE VALUE FOR A FLAG "IFLAGR" TO DETERMINE IF RNEW IS TO BE 00000400 C CALCULATED (1) OR NOT (0) AND A FLAG "IFLAGS" TO DETERMINE IF 00000410 C SMOOTHED ESTIMATES FOR INTERMEDIATE X VALUES ARE TO BE CALCULATED 00000420 C (1) OR NOT (0). 00000430 C 00000440 CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00000450 CALL SHIFT(N,N,QZERO,QO) 00000460 CALL SHIFT(N,1,PZERO,PO) 00000470 CALL SHIFT(1,1,RZERO,RO) 00000480 C 00000490 C ENTERING THE MAIN DO LOOP FOR GENERATING Q,P,AND R FOR 00000500 C SUCCESSIVE TIMES T = 1,TCAP USING EQS.(24), (26), AND (27). 00000510 C 00000520 DO 50 T=1,TCAP 00000530 CALL MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00000540 DO 5 I=1,MOBS 00000550 YY(I,T) = Y(I,1) 00000560 5 CONTINUE 00000570 C 00000580 C GETTING U=HT*M*H + QO IN EQ.(28) 00000590 C 00000600 CALL MUL(MOBS,MOBS,N,M,H,AA) 00000610 CALL TRANS(MOBS,N,H,HT) 00000620 CALL MUL(N,MOBS,N,HT,AA,BB) 00000630 CALL ADD(N,N,BB,QO,U) 00000640 C 00000650 C GETTING C=FT*D 00000660 C 00000670 CALL TRANS(N,N,F,AA) 00000680 CALL MUL(N,N,N,AA,D,C) 00000690 C 00000700 C GETTING W=AMU*C*F+U 00000710 C 00000720 CALL MUL(N,N,N,C,F,AA) 00000730 CALL MULCON(N,N,AMU,AA,BB) 00000740 CALL ADD(N,N,BB,U,W) 00000750 C 00000760 C GETTING V=WINV IN EQ.(17) 00000770 C 00000780 CALL INV(N,W,V) 00000790 C 00000800 C GETTING E = (Y-B) 00000810 C 00000820 CALL SUB(MOBS,1,Y,B,E) 00000830 C 00000840 C GETTING Z = HT*M*E + PO IN EQ.(29) 00000850 C 00000860 CALL MUL(MOBS,MOBS,1,M,E,AA) 00000870 CALL MUL(N,MOBS,1,HT,AA,BB) 00000880 CALL ADD(N,1,BB,PO,Z) 00000890 C 00000900 C GETTING G = AMU*V*C IN EQ.(20) 00000910 C 00000920 CALL MUL(N,N,N,V,C,AA) 00000930 CALL MULCON(N,N,AMU,AA,G) 00000940 IF(IFLAGS.EQ.0) GO TO 110 00000950 C 00000960 C STORE G FOR CALCULATION OF SMOOTHED ESTIMATES 00000970 C 00000980 DO 10 I=1,N 00000990 DO 20 J=1,N 00001000 GG(I,J,T)=G(I,J) 00001010 20 CONTINUE 00001020 10 CONTINUE 00001030 110 CONTINUE 00001040 C 00001050 C GETTING QNEW = AMU*D*(I-F*G) IN EQ.(24) 00001060 C 00001070 CALL MUL(N,N,N,F,G,AA) 00001080 CALL IDEN(N,BB) 00001090 CALL SUB(N,N,BB,AA,CC) 00001100 CALL MUL(N,N,N,D,CC,DD) 00001110 CALL MULCON(N,N,AMU,DD,QNEW) 00001120 C 00001130 C GETTING PNEW = GT*Z+QNEWT*A IN EQ.(26) 00001140 C 00001150 CALL TRANS(N,N,G,AA) 00001160 CALL MUL(N,N,1,AA,Z,BB) 00001170 CALL TRANS(N,N,QNEW,CC) 00001180 CALL MUL(N,N,1,CC,A,DD) 00001190 CALL ADD(N,1,BB,DD,PNEW) 00001200 C 00001210 C GETTING S = V*(Z - AMU*C*A) IN EQ.(19) 00001220 C 00001230 CALL MUL(N,N,1,C,A,BB) 00001240 CALL MULCON(N,1,AMU,BB,CC) 00001250 CALL SUB(N,1,Z,CC,DD) 00001260 CALL MUL(N,N,1,V,DD,S) 00001270 IF(IFLAGS.EQ.0) GO TO 210 00001280 C 00001290 C STORE S FOR CALCULATION OF SMOOTHED ESTIMATES 00001300 C 00001310 DO 30 I=1,N 00001320 SS(I,T)=S(I,1) 00001330 30 CONTINUE 00001340 210 CONTINUE 00001350 IF(IFLAGR.EQ.0) GO TO 310 00001360 C 00001370 C GETTING RNEW = RO + ET*M*E + AMU*AT*D*A - ST*W*S IN EQ.(27) 00001380 C 00001390 CALL MUL(MOBS,MOBS,1,M,E,AA) 00001400 CALL TRANS(MOBS,1,E,BB) 00001410 CALL MUL(1,MOBS,1,BB,AA,CC) 00001420 CALL ADD(1,1,RO,CC,DD) 00001430 CALL MUL(N,N,1,D,A,EE) 00001440 CALL TRANS(N,1,A,FF) 00001450 CALL MUL(1,N,1,FF,EE,HH) 00001460 CALL MULCON(1,1,AMU,HH,OO) 00001470 CALL ADD(1,1,DD,OO,PP) 00001480 CALL MUL(N,N,1,W,S,QQ) 00001490 CALL TRANS(N,1,S,RR) 00001500 CALL MUL(1,N,1,RR,QQ,TT) 00001510 CALL SUB(1,1,PP,TT,RNEW) 00001520 310 CONTINUE 00001530 IF(T.EQ.TCAP) GO TO 50 00001540 C 00001550 C UPDATING QO,PO, AND RO 00001560 C 00001570 CALL SHIFT(N,N,QNEW,QO) 00001580 CALL SHIFT(N,1,PNEW,PO) 00001590 IF(IFLAGR.EQ.0) GO TO 50 00001600 CALL SHIFT(1,1,RNEW,RO) 00001610 50 CONTINUE 00001620 C 00001630 C GETTING THE FLS FILTER ESTIMATE FOR XTCAP = UINV*Z IN EQ.(30) 00001640 C 00001650 CALL INV(N,U,AA) 00001660 CALL MUL(N,N,1,AA,Z,XTCAP) 00001670 DO 65 I=1,N 00001680 X(I,TCAP)=XTCAP(I,1) 00001690 65 CONTINUE 00001700 IF (IFLAGS.EQ.1) GOTO 410 00001710 C 00001720 C PRINTING OUT THE FLS FILTER ESTIMATE FOR XTCAP 00001730 C 00001740 CALL OUTPUT(TCAP,N,X,TRUEX) 00001750 IF(IFLAGS.EQ.0) GOTO 510 00001760 410 CONTINUE 00001770 C 00001780 C GETTING SMOOTHED ESTIMATES FOR X1,..., XTCAP-1 IN EQS.(33A) 00001790 C 00001800 TCAP1=TCAP-1 00001810 DO 70 T=1,TCAP1 00001820 L=TCAP-T 00001830 DO 80 I=1,N 00001840 X(I,L)=SS(I,L) 00001850 DO 90 J=1,N 00001860 X(I,L)=X(I,L)+GG(I,J,L)*X(J,L+1) 00001870 90 CONTINUE 00001880 80 CONTINUE 00001890 70 CONTINUE 00001900 C 00001910 C PRINTING OUT THE FLS ESTIMATES FOR X1,...,XTCAP 00001920 C 00001930 DO 150 T=1,TCAP 00001940 CALL OUTPUT(T,N,X,TRUEX) 00001950 150 CONTINUE 00001960 C VALIDATION TEST: HOW WELL DO THE FLS ESTIMATES SATISFY THE 00001970 C FIRST-ORDER CONDITIONS FOR THE COST MINIMIZATION PROBLEM (5) 00001980 CALL FOCTST(X,YY) 00001990 510 CONTINUE 00002000 STOP 00002010 END 00002020 C 00002030 C MATRIX SUBROUTINES FOR ADDITION, MULTIPLICATION, TRANSPOSITION, 00002040 C SUBTRACTION, INVERSION, MULTIPLICATION BY A SCALAR, SHIFT, AND 00002050 C FORMATION OF AN IDENTITY MATRIX 00002060 C 00002070 C OBTAINING THE SUM C=A+B OF TWO NROW X MCOL MATRICES A AND B 00002080 C 00002090 SUBROUTINE ADD(NROW,MCOL,A,B,C) 00002100 IMPLICIT REAL*8(A-H,O-Z) 00002110 DIMENSION A(15,15),B(15,15),C(15,15) 00002120 DO 10 I=1,NROW 00002130 DO 20 J=1,MCOL 00002140 C(I,J)=A(I,J)+B(I,J) 00002150 20 CONTINUE 00002160 10 CONTINUE 00002170 RETURN 00002180 END 00002190 C 00002200 C OBTAINING THE PRODUCT C=A*B OF AN NROW X L MATRIX A AND AN 00002210 C L X MCOL MATRIX B 00002220 C 00002230 SUBROUTINE MUL(NROW,L,MCOL,A,B,C) 00002240 IMPLICIT REAL*8(A-H,O-Z) 00002250 DIMENSION A(15,15),B(15,15),C(15,15) 00002260 DO 10 I=1,NROW 00002270 DO 20 J=1,MCOL 00002280 SUM=0.0D+00 00002290 DO 30 K=1,L 00002300 SUM=SUM+A(I,K)*B(K,J) 00002310 30 CONTINUE 00002320 C(I,J)=SUM 00002330 20 CONTINUE 00002340 10 CONTINUE 00002350 RETURN 00002360 END 00002370 C 00002380 C OBTAINING THE TRANSPOSE B OF AN NROW X MCOL MATRIX A 00002390 C 00002400 SUBROUTINE TRANS(NROW,MCOL,A,B) 00002410 IMPLICIT REAL*8(A-H,O-Z) 00002420 DIMENSION A(15,15),B(15,15) 00002430 DO 10 I=1,NROW 00002440 DO 20 J=1,MCOL 00002450 B(J,I)=A(I,J) 00002460 20 CONTINUE 00002470 10 CONTINUE 00002480 RETURN 00002490 END 00002500 C 00002510 C OBTAINING THE DIFFERENCE C=A-B BETWEEN NROW X MCOL MATRICES 00002520 C A AND B 00002530 C 00002540 SUBROUTINE SUB(NROW,MCOL,A,B,C) 00002550 IMPLICIT REAL*8(A-H,O-Z) 00002560 DIMENSION A(15,15),B(15,15),C(15,15) 00002570 DO 10 I=1,NROW 00002580 DO 20 J=1,MCOL 00002590 C(I,J)=A(I,J)-B(I,J) 00002600 20 CONTINUE 00002610 10 CONTINUE 00002620 RETURN 00002630 END 00002640 C 00002650 C OBTAINING THE INVERSE C OF A K X K MATRIX A 00002660 C 00002670 SUBROUTINE INV(K,A,C) 00002680 IMPLICIT REAL*8(A-H,O-Z) 00002690 DIMENSION A(15,15),B(15,30),C(15,15) 00002700 DO 5 J=1,K 00002710 DO 6 I=1,K 00002720 B(I,J)=A(I,J) 00002730 6 CONTINUE 00002740 5 CONTINUE 00002750 K2=K*2 00002760 DO 7 J=1,K 00002770 DO 8 I=1,K 00002780 B(I,K+J)=0.0D+00 00002790 IF(I.EQ.J) B(I,K+J)=1.0D+00 00002800 8 CONTINUE 00002810 7 CONTINUE 00002820 C 00002830 C THE PIVOT OPERATION STARTS HERE 00002840 C 00002850 DO 9 L=1,K 00002860 PIVOT = B(L,L) 00002870 DO 13 J=L,K2 00002880 B(L,J)=B(L,J)/PIVOT 00002890 13 CONTINUE 00002900 C 00002910 C TO IMPROVE THE ROWS 00002920 C 00002930 DO 14 I=1,K 00002940 IF(I.EQ.L) GO TO 14 00002950 AIL=B(I,L) 00002960 DO 15 J=L,K2 00002970 B(I,J)=B(I,J)-AIL*B(L,J) 00002980 15 CONTINUE 00002990 14 CONTINUE 00003000 9 CONTINUE 00003010 DO 45 I=1,K 00003020 DO 46 J=1,K 00003030 C(I,J)=B(I,K+J) 00003040 46 CONTINUE 00003050 45 CONTINUE 00003060 RETURN 00003070 END 00003080 C 00003090 C OBTAINING THE PRODUCT C*A OF A SCALAR C AND AN NROW X MCOL 00003100 C MATRIX A 00003110 C 00003120 SUBROUTINE MULCON(NROW,MCOL,C,A,CA) 00003130 IMPLICIT REAL*8(A-H,O-Z) 00003140 DIMENSION A(15,15),CA(15,15) 00003150 DO 10 I=1,NROW 00003160 DO 20 J=1,MCOL 00003170 CA(I,J)=C*A(I,J) 00003180 20 CONTINUE 00003190 10 CONTINUE 00003200 RETURN 00003210 END 00003220 C 00003230 C PUTTING AN NROW X MCOL MATRIX A INTO AN NROW X MCOL MATRIX B 00003240 C 00003250 SUBROUTINE SHIFT(NROW,MCOL,A,B) 00003260 IMPLICIT REAL*8(A-H,O-Z) 00003270 DIMENSION A(15,15),B(15,15) 00003280 DO 10 I=1,NROW 00003290 DO 20 J=1,MCOL 00003300 B(I,J)=A(I,J) 00003310 20 CONTINUE 00003320 10 CONTINUE 00003330 RETURN 00003340 END 00003350 C 00003360 C FORMING THE N X N IDENTITY MATRIX E 00003370 C 00003380 SUBROUTINE IDEN(N,E) 00003390 IMPLICIT REAL*8(A-H,O-Z) 00003400 DIMENSION E(15,15) 00003410 ZERO=0.0D+00 00003420 ONE=1.0D+00 00003430 DO 10 I=1,N 00003440 DO 20 J=1,N 00003450 E(I,J)=ZERO 00003460 20 CONTINUE 00003470 10 CONTINUE 00003480 DO 30 L=1,N 00003490 E(L,L)=ONE 00003500 30 CONTINUE 00003510 RETURN 00003520 END 00003530 C 00003540 SUBROUTINE INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00003550 IMPLICIT REAL*8(A-H,O-Z) 00003560 INTEGER TCAP 00003570 DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00003580 AMU = 1.0D+00 00003590 TCAP = 30 00003600 N = 2 00003610 MOBS = 1 00003620 DO 10 J = 1,N 00003630 DO 20 I = 1,N 00003640 QZERO(I,J) = 0.0D+00 00003650 PZERO(I,J) = 0.0D+00 00003660 RZERO(I,J) = 0.0D+00 00003670 20 CONTINUE 00003680 10 CONTINUE 00003690 IFLAGR=1 00003700 IFLAGS=1 00003710 RETURN 00003720 END 00003730 C 00003740 C 00003750 SUBROUTINE MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00003760 IMPLICIT REAL*8(A-H,O-Z) 00003770 REAL*8 M 00003780 INTEGER T,TCAP 00003790 DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15),M(15,15) 00003800 DIMENSION Y(15,15),TRUEX(15,110),ZERO(15,15) 00003810 DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00003820 C 00003830 C REGIME SHIFT EXPERIMENT WITH POSSIBLY NOISY OBSERVATIONS 00003840 C 00003850 CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00003860 SIGMA = 0.00D+00 00003870 DO 10 I=1,15 00003880 DO 20 J=1,15 00003890 ZERO(I,J) = 0.0D+00 00003900 20 CONTINUE 00003910 10 CONTINUE 00003920 CALL IDEN(N,F) 00003930 CALL SHIFT(N,1,ZERO,A) 00003940 H(1,1)=1.0D+00 00003950 H(1,2)=1.0D+00 00003960 AT=DFLOAT(T) 00003970 IF(T.EQ.1) GO TO 200 00003980 H(1,1)=DSIN(10.0D+00+(AT))+0.01D+00 00003990 H(1,2)=DCOS(10.0D+00+(AT)) 00004000 200 CONTINUE 00004010 CALL SHIFT(MOBS,1,ZERO,B) 00004020 CALL IDEN(N,D) 00004030 CALL IDEN(MOBS,M) 00004040 IF (T.GT.15) GOTO 150 00004050 TRUEX(1,T) = 2.0D+00 00004060 TRUEX(2,T) = 3.0D+00 00004070 GOTO 175 00004080 150 TRUEX(1,T) = 4.0D+00 00004090 TRUEX(2,T) = 5.0D+00 00004100 175 CONTINUE 00004110 C THE DISCREPANCY TERM UU IN THE NEXT LINE IS SET TO ZERO; C IT CAN BE REPLACED BY A CALL TO A PSEUDO-RANDOM NUMBER GENERATOR C TO TEST GFLS IN THE PRESENCE OF NOISE IN THE OBSERVATIONS UU = 0.0D+00 00004120 Y(1,1)=H(1,1)*TRUEX(1,T) + H(1,2)*TRUEX(2,T) + UU 00004130 RETURN 00004140 END 00004150 C 00004160 SUBROUTINE OUTPUT(T,N,X,TRUEX) 00004170 IMPLICIT REAL*8(A-H,O-Z) 00004180 INTEGER T 00004190 DIMENSION X(15,110),TRUEX(15,110) 00004200 L = T 00004210 WRITE(6,100) L,(X(I,L),I=1,N) 00004220 100 FORMAT(1H0,'TIME EQUALS',I3/1X,'FLS ESTIMATES',7X,2D25.10) 00004230 WRITE(6,200) (TRUEX(I,L),I=1,N) 00004240 200 FORMAT(1X,'TRUE X VALUES',7X,2D25.10) 00004250 RETURN 00004260 END 00004270 C 00004280 C VALIDATION TEST: HOW WELL DO THE FLS ESTIMATES SATISFY THE 00004290 C FIRST-ORDER CONDITIONS FOR THE COST MINIMIZATION PROBLEM (5) 00004300 C 00004310 SUBROUTINE FOCTST(X,YY) 00004320 IMPLICIT REAL*8(A-H,O-Z) 00004330 INTEGER T,TP1,TCAP,TCAP1 00004340 REAL*8 M,MH 00004350 DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00004360 DIMENSION XT(15,15),X(15,110),XTT(15,15),E(15,15) 00004370 DIMENSION PZEROT(15,15),EE(15,15),CO(15,15),YT(15,15),YY(15,110) 00004380 DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15),M(15,15) 00004390 DIMENSION Y(15,15),TRUEX(15,110) 00004400 DIMENSION MH(15,15),EM(15,15),EMT(15,15),W(15,15),XTP1(15,15) 00004410 DIMENSION ED(15,15),EDT(15,15),U(15,15),V(15,15),FOCD(15,15) 00004420 C = -1.0D+00 00004430 C FORM THE STATE VECTOR FOR TIME T = 1 00004440 CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00004450 DO 100 I=1,N 00004460 XT(I,1) = X(I,1) 00004470 100 CONTINUE 00004480 C FORM THE INITIAL INCREMENTAL COST CO = -(X1'QO - PO') 00004490 CALL TRANS(N,1,XT,XTT) 00004500 CALL MUL(1,N,N,XTT,QZERO,E) 00004510 CALL TRANS(N,1,PZERO,PZEROT) 00004520 CALL SUB(1,N,E,PZEROT,EE) 00004530 CALL MULCON(1,N,C,EE,CO) 00004540 C DO LOOP FOR THE SEQUENTIAL CHECK OF THE FOC FOR T=1,TCAP 00004550 DO 200 T=1,TCAP 00004560 C FORM THE TIME-T STATE VECTOR XT 00004570 DO 300 I=1,N 00004580 XT(I,1) = X(I,T) 00004590 300 CONTINUE 00004600 C FORM THE TIME-T OBSERVATION VECTOR YT 00004610 DO 400 J=1,MOBS 00004620 YT(J,1) = YY(J,T) 00004630 400 CONTINUE 00004640 CALL MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00004650 C FORM W = (YT - H(T)XT - B(T))'M(T)H(T) 00004660 CALL MUL(MOBS,MOBS,N,M,H,MH) 00004670 CALL RME(N,MOBS,YT,XT,H,B,EM) 00004680 CALL TRANS(MOBS,1,EM,EMT) 00004690 CALL MUL(1,MOBS,N,EMT,MH,W) 00004700 IF(T.EQ.TCAP) GOTO 600 00004710 C FORM THE TIME-T+1 STATE VECTOR XTP1 00004720 TP1 = T + 1 00004730 DO 500 I=1,N 00004740 XTP1(I,1) = X(I,TP1) 00004750 500 CONTINUE 00004760 C FORM U = AMU*(XTP1 - F(T)XT - A(T))'*D(T) 00004770 CALL RDE(N,XTP1,XT,F,A,ED) 00004780 CALL TRANS(N,1,ED,EDT) 00004790 CALL MUL(1,N,N,EDT,D,E) 00004800 CALL MULCON(1,N,AMU,E,U) 00004810 C FORM V = U*F 00004820 CALL MUL(1,N,N,U,F,V) 00004830 GOTO 800 00004840 600 CONTINUE 00004850 DO 700 I=1,N 00004860 V(1,I) = 0.0D+00 00004870 700 CONTINUE 00004880 800 CONTINUE 00004890 C DETERMINE THE FOC DISCREPANCIES FOR TIME T 00004900 C GIVEN BY FOCD = CO + V + W 00004910 CALL ADD(1,N,CO,V,E) 00004920 CALL ADD(1,N,E,W,FOCD) 00004930 C PRINT OUT THE FOC DISCREPANCIES FOCD FOR TIME T 00004940 WRITE (6,36) T 00004950 36 FORMAT(1H0,'FOC DISCREPANCIES FOR TIME',I3) 00004960 WRITE (6,37) (FOCD(1,I),I=1,N) 00004970 37 FORMAT(1X,13D10.2) 00004980 C UPDATE THE INITIAL INCREMENTAL COST CO 00004990 CALL MULCON(1,N,C,U,CO) 00005000 200 CONTINUE 00005010 RETURN 00005020 END 00005030 C 00005040 C SUBROUTINE FOR EVALUATING THE MEASUREMENT SPECIFICATION ERROR 00005050 C EM = (YT - H(T)XT - B(T)) FOR TIME T 00005060 C 00005070 SUBROUTINE RME(N,MOBS,YT,XT,H,B,EM) 00005080 IMPLICIT REAL*8(A-H,O-Z) 00005090 DIMENSION YT(15,15),XT(15,15),H(15,15),B(15,15),EM(15,15) 00005100 DIMENSION HX(15,15),HXPB(15,15) 00005110 CALL MUL(MOBS,N,1,H,XT,HX) 00005120 CALL ADD(MOBS,1,HX,B,HXPB) 00005130 CALL SUB(MOBS,1,YT,HXPB,EM) 00005140 RETURN 00005150 END 00005160 C 00005170 C SUBROUTINE FOR EVALUATING THE DYNAMIC SPECIFICATION ERROR 00005180 C ED = (XTP1 - F(T)XT - A(T)) FOR TIME T 00005190 C 00005200 SUBROUTINE RDE(N,XTP1,XT,F,A,ED) 00005210 IMPLICIT REAL*8(A-H,O-Z) 00005220 DIMENSION XTP1(15,15),XT(15,15),F(15,15),A(15,15),ED(15,15) 00005230 DIMENSION FXT(15,15),FXTPA(15,15) 00005240 CALL MUL(N,N,1,F,XT,FXT) 00005250 CALL ADD(N,1,FXT,A,FXTPA) 00005260 CALL SUB(N,1,XTP1,FXTPA,ED) 00005270 RETURN 00005280 END 00005290