Квазирекуррентный алгоритм идентификации .

Программа идентификации по методу квадратного корня.

Программа №1 /До 8 коэффициентов а, и 8 коэффициентов b, и 8 шагов max запаздывания, т.е. max сдвиг аргумента до 16 /

DIMENSION Y(8),U(16),A(8),B(8),GM(153)

DIMENSION YP(100),UP(100)

READ(5,5) L,NDMAX /Число пар;вход-выход-L; max запаздывание: NDMAX/

5 FORMAT(I3,I2)

WRITE(6,25) L,NDMAX

25 FORMAT(1H /10X,'Исходные данные для идентификации'//

*10X,'L=',I2,' NDMAX=',I2/1H )

READ(5,15) (UP(I),YP(I),I=1,L) /UP-вход, YP-выход, I-дискретное время до: L/

15 FORMAT(2F8.4)

WRITE(6,35) (UP(I),YP(I),I=1,L)

35 FORMAT(1H /12X,'BXOД',9X,'BЫXOД'//50(10X,F8.4,5X,F8.4/))

NDMAX=NDMAX+1 /Сдвиг на единицу, чтобы вариант нулевого сдвига был при индексе:1 /

DO 70 K=1,NDMAX /Цикл многократного решения для разных запаздываний/

ND=K-1

CALL MSF(1,3,ND,YN,Y,U,A,B,CGL,GM) /Начальное заполнение массивов/

DO 10 I=1,3

A(I)=0.0

10 B(I)=0.0

CGL=0.0

CKO=0.0

DO 60 M=1,2 /М=1 – поиск параметров, М=2 – оценка найденной модели /

DO 60 I=1,L /Главный цикл считывания пар: вход-выход/

YN=YP(I)

IF(I.LT.(4+ND)) GOTO 30 /Набор первого варианта. j здесь для 3-х ‘а’ и 3-х ‘b’. Т.е. надо 3 пары вход-выход чтобы начать итерации/

IF(M.EQ.2) GOTO 20

CALL MSF(2,3,ND,YN,Y,U,A,B,CGL,GM) /Рабочий шаг идентификации/

GOTO 30

20 SKO=YN-B(1)*U(ND+1)-B(2)*U(ND+2)-B(3)*U(ND+3)

*+A(1)*Y(1)+A(2)*Y(2)+A(3)*Y(3)-CGL

CKO=CKO+SKO*SKO

30 DO 40 J=1,2

40 Y(4-J)=Y(3-J)

DO 50 J=1,15

50 U(17-J)=U(16-J)

Y(1)=YN

60 U(1)=UP(I)

CKO=SQRT(CKO)/(L-ND-3)

WRITE(6,45) ND,(I,A(I),I=1,3),(I,B(I),I=1,3),CGL,CKO

45 FORMAT(1H /10X,'PЕЗУЛЬТАТ ИДЕНТИФИКАЦИИ ND=',I2//

*10X,3(' A',I1,' = ',G12.5)/10X,3(' B',I1,' = ',G12.5)/10X,

*' C0= ',G12.5//10X,' CKO=',G12.5/1H )

70 CONTINUE

STOP

END

SUBROUTINE MSF(INIT,M,ND,YN,Y,U,A,B,CGL,GM)

DIMENSION Y(8),U(16),A(8),B(8),GM(153)

DIMENSION TETA(17),PSI(17),GV(17)

INU(I,J)=I*(I-1)/2+J /Формула одномерного индекса в треугольном массиве: i,j/

N=2*M+1

IF(INIT.NE.1) GOTO 20

10 DO 14 I=1,N

DO 12 J=1,I

12 GM(INU(I,J))=0.0

14 GM(INU(I,I))=1.0E+5

RETURN

20 DO 25 I=1,M

TETA(2*I-1)=A(I)

TETA(2*I)=B(I)

PSI(2*I-1)=-Y(I)

25 PSI(2*I)=U(I+ND)

TETA(N)=CGL

PSI(N)=1.0

SIG=1.0

DO 40 L=1,N

HL=0.0

DO 32 I=1,L

32 HL=HL+GM(INU(L,I))*PSI(I)

SIG1=SIG

SIG=SIG+HL*HL

PHISS=SQRT(SIG1/SIG)

GV(L)=HL*GM(INU(L,L))

GM(INU(L,L))=PHISS*GM(INU(L,L))

IF(L-1) 40,40,35

35 LUA=L-1

DO 39 J=1,LUA

GMLJ=GM(INU(L,J))

GM(INU(L,J))=PHISS*(GMLJ-HL*GV(J)/SIG1)

39 GV(J)=GV(J)+HL*GMLJ

40 CONTINUE

EPSN=YN

DO 50 I=1,N

50 EPSN=EPSN-PSI(I)*TETA(I)

EPSN=EPSN/SIG

DO 52 I=1,N

52 TETA(I)=TETA(I)+EPSN*GV(I)

60 DO 62 I=1,M

A(I)=TETA(2*I-1)

62 B(I)=TETA(2*I)

CGL=TETA(N)

RETURN

END

Исходные данные для идентификации

L= 60 NDMAX= 8

BXOД BЫXOД

-1.4710 5.1980

-.5785 6.3200

1.0460 5.8260

-.5766 3.3380

-.4351 -3.6970

-1.2630 -8.3750

.7064 -5.0110

.9887 4.0708

.0995 -1.7370

.8031 -2.6910

.3668 -6.9860

-1.1940 1.7610

.1346 5.5190

.6195 1.9070

.0273 4.4190

-.1278 2.9480

-.3297 -5.2780

-.8129 -.7508

.6869 2.9880

1.1880 .9287

-1.2950 -.4464

.5024 -1.7820

-.6517 -4.5210

-.2448 2.2890

-.7592 6.6120

.5688 -4.7900

-.9917 1.1560

-1.4290 -2.8820

.2751 -1.9900

-.6137 -4.2660

-.3316 1.7690

-.3619 -4.4260

-.1892 -8.3220

-.0478 -.7158

.3277 -3.1140

.9536 -2.4510

.8451 -2.3940

.2682 -1.5260

.9878 -.5946

1.5210 1.5110

.3804 5.1710

.3346 5.5450

-1.3520 2.6960

.3365 5.5450

-1.5840 8.9980

.1512 4.1480

.3082 2.5990

.4029 -6.1560

-.0488 .0363

-.2819 -7.8120

-2.0330 -1.2730

.7191 1.3400

.9858 2.3880

.0341 .3558

.6283 -1.3570

-.7133 -10.5200

.5209 .8815

.1005 5.3340

-.7103 1.5490

.1024 3.4550

Результат идентификации ND= 0

A1 = -.32046 А2 = .53016E-01 A3 = .66496E-01

B1 = -.73783E-01 B2 = -.18001 B3 = -.22090E-01

C0 = -.20612

CKO= .52723

Результат идентификации ND= 1

A1 = -.30979 A2 = .60306E-01 A3 = .73355E-01

B1 = -.12810 B2 = .74783E-01 B3 = .21790

C0= -.24128

CKO= .53447

Результат идентификации ND= 2

A1 = -.26002 A2 = .16607E-01 A3 = .55256E-03

B1 = -.21285E-05 B2 = .16749E-04 B3 = 4.9996

C0= .20490E-03

CKO= .12523E-03

Результат идентификации ND= 3

A1 = -.31070 A2 = .29786E-01 A3 = -.30287E-03

B1 = .12813E-04 B2 = 4.9995 B3 = -.25341

C0= .20458E-03

CKO= .12525E-03

Результат идентификации ND= 4

A1 = -.15616 A2 = -.15956E-01 A3 = .37029E-02

B1 = 4.9995 B2 = .51919 B3 = -.27818E-01

C0= .24978E-03

CKO= .12505E-03

Результат идентификации ND= 5

A1 = -114.11 A2 = 42.109 A3 = -4.5782

B1 = -569.08 B2 = 62.597 B3 = 2.8433

C0= -.90414E-01

CKO= .54420

Результат идентификации ND= 6

A1 = -.30302 A2 = 138.37 A3 = -23.096

B1 = 691.75 B2 = 64.289 B3 = 5.1893

C0 = .24715E-01

CKO= .55133

Результат идентификации ND= 7

A1 = -.29190 A2 = -.38948E-02 A3 = -15.798

B1 = -79.128 B2 = -20.721 B3 = -4.1379

C0= .47414E-01

CKO= .55604

Результат идентификации ND= 8

A1 = -.27850 A2 = .16188E-01 A3 = -.24854E-01

B1 = -.30679 B2 = -.43860 B3 = .66448

C0= .16573

CKO= .55087

Stop - Program terminated.

Программа №2

DIMENSION Y(8),U(16),A2(8),A(10,8),B2(8),B(10,8),GM(153)

DIMENSION YP(100),UP(100),CGL(10),CKO(10)

READ(5,5) L1,NDMAX

5 FORMAT(I3,I2)

WRITE(6,25) L1,NDMAX

25 FORMAT(1H /10X,'ИCXOДHЫE ДAHHЫE ДЛЯ ИДEHTИФИKAЦИИ'//

*10X,'L1=',I2,' NDMAX=',I2/1H )

READ(5,15) (UP(I),YP(I),I=1,L1)

15 FORMAT(2F8.4)

OPEN(3,FILE = 'LAB.TXT')

WRITE(6,35) (UP(I),YP(I),I=1,L1)

35 FORMAT(1H /12X,'BXOД',9X,'BЫXOД'//50(10X,F8.4,5X,F8.4/))

NDMAX=NDMAX+1

DO 70 K=1,NDMAX

ND=K-1

CALL MSF(1,3,ND,YN,Y,U,A,B,CGL,GM)

DO 10 I=1,3

A(K,I)=0.0

10 B(K,I)=0.0

CGL(K)=0.0

CKO(K)=0.0

DO 60 M=1,2

DO 60 I=1,L1

YN=YP(I)

IF(I.LT.(4+ND)) GOTO 30

IF(M.EQ.2) GOTO 20

CALL MSF(2,3,ND,YN,Y,U,A,B,CGL,GM)

GOTO 30

20 SKO=YN-B(K,1)*U(ND+1)-B(K,2)*U(ND+2)-B(K,3)*U(ND+3)

*+A(K,1)*Y(1)+A(K,2)*Y(2)+A(K,3)*Y(3)-CGL(K)

CKO(K)=CKO(K)+SKO*SKO

30 DO 40 J=1,2

40 Y(4-J)=Y(3-J)

DO 50 J=1,15

50 U(17-J)=U(16-J)

Y(1)=YN

60 U(1)=UP(I)

CKO(K)=SQRT(CKO(K))/(L1-ND-3)

WRITE (6,135) ND,(I,A(K,I),I=1,3),(I,B(K,I),I=1,3),CGL(K),CKO(K)

135 FORMAT (/10X,'rezultat of identification ND=',I2//

*10X,3(' A',I1,'=',G12.5)/10X,3(' B',I1,'=',G12.5)/10X,

*' C0= ',G12.5//10X,' CKO=',G12.5/)

70 CONTINUE

CKOMIN=CKO(1)

DO 80 I=1,NDMAX

IF (CKO(I).GE.CKOMIN) GOTO 80

CKOMIN=CKO(I)

IMIN=I

80 CONTINUE

WRITE (6,140) IMIN-1,CKOMIN

140 FORMAT (/10X,'CKOmin=CKO(',I1,')=',G12.5)

WRITE (6,145) IMIN-1,(I,A(IMIN,I),I=1,3),(I,B(IMIN,I),I=1,3),

*CGL(IMIN),CKO(IMIN)

145 FORMAT (/15X,'Optimal model: ND=',I2,//10X,3(' A',I1,'=',G12.5)

*/10X,3(' B',I1,'=',G12.5)/10X,' C0= ',G12.5//10X,' CKO=',G12.5/)

L=60

WRITE (6,160) L1+1

160 FORMAT (/10X,'ishodn. data for identification'//

*10X,'L2=',I2/)

READ (3,165) (UP(I),I=L1+1,L)

165 FORMAT (F7.4)

DO 90 I=1,3

A2(I)=A(IMIN,I)

90 B2(I)=B(IMIN,I)

CGL2=CGL(IMIN)

ND=IMIN-1

DO 95 K=(L1+1),L

95 YP(K)=-A2(1)*YP(K-1)-A2(2)*YP(K-2)-A2(3)*YP(K-3)+B2(1)*

*UP(K-ND-1)+B2(2)*UP(K-ND-2)+B2(3)*UP(K-ND-3)+CGL2

WRITE (6,170)

170 FORMAT (/8X,'Output for model optimal',/9X,'INPUT',8X,'OUTPUT',/)

WRITE (6,175) (UP(I),YP(I),I=1,L)

175 FORMAT (5X,F9.4,5X,F9.4)

STOP

END

SUBROUTINE MSF(INIT,M,ND,YN,Y,U,A,B,CGL,GM)

DIMENSION Y(8),U(16),A(10,8),B(10,8),GM(153),CGL(10)

DIMENSION TETA(17),PSI(17),GV(17)

INU(I,J)=I*(I-1)/2+J

N=2*M+1

IF (INIT.NE.1) GOTO 20

10 DO 14 I=1,N

DO 12 J=1,I

12 GM(INU(I,J))=0.0

14 GM(INU(I,I))=1.0E+5

RETURN

20 DO 25 I=1,M

TETA(2*I-1)=A(ND+1,I)

TETA(2*I)=B(ND+1,I)

PSI(2*I-1)=-Y(I)

25 PSI(2*I)=U(I+ND)

TETA(N)=CGL(ND+1)

PSI(N)=1.0

SIG=1.0

DO 40 L=1,N

HL=0.0

DO 32 I=1,L

32 HL=HL+GM(INU(L,I))*PSI(I)

SIG1=SIG

SIG=SIG+HL*HL

PHISS=SQRT(SIG1/SIG)

GV(L)=HL*GM(INU(L,L))

GM(INU(L,L))=PHISS*GM(INU(L,L))

IF (L-1) 40,40,35

35 LUA=L-1

DO 39 J=1,LUA

GMLJ=GM(INU(L,J))

GM(INU(L,J))=PHISS*(GMLJ-HL*GV(J)/SIG1)

39 GV(J)=GV(J)+HL*GMLJ

40 CONTINUE

EPSN=YN

DO 50 I=1,N

50 EPSN=EPSN-PSI(I)*TETA(I)

EPSN=EPSN/SIG

DO 52 I=1,N

52 TETA(I)=TETA(I)+EPSN*GV(I)

60 DO 62 I=1,M

A(ND+1,I)=TETA(2*I-1)

62 B(ND+1,I)=TETA(2*I)

CGL(ND+1)=TETA(N)

RETURN

END

Исходные данные для идентификации

L1=30 NDMAX= 8

BXOД BЫXOД

-1.4710 5.1980

-.5785 6.3200

1.0460 5.8260

-.5766 3.3380

-.4351 -3.6970

-1.2630 -8.3750

.7064 -5.0110

.9887 4.0708

.0995 -1.7370

.8031 -2.6910

.3668 -6.9860

-1.1940 1.7610

.1346 5.5190

.6195 1.9070

.0273 4.4190

-.1278 2.9480

-.3297 -5.2780

-.8129 -.7508

.6869 2.9880

1.1880 .9287

-1.2950 -.4464

.5024 -1.7820

-.6517 -4.5210

-.2448 2.2890

-.7592 6.6120

.5688 -4.7900

-.9917 1.1560

-1.4290 -2.8820

.2751 -1.9900

-.6137 -4.2660

rezultat of identification ND= 0

A1= -.16128 A2= .23082 A3= .12149

B1= .42131 B2= .60180E-01 B3= -.99698

C0= -.61669

CKO= .68526

rezultat of identification ND= 1

A1= -.39983E-01 A2= .27186 A3= .20675

B1= .90560E-01 B2= -.92743 B3= -1.0363

C0= -.92056

CKO= .68006

rezultat of identification ND= 2

A1= -.25998 A2= .16633E-01 A3= .58827E-03

B1= .12900E-04 B2= -.15851E-03 B3= 4.9999

C0= .23152E-03

CKO= .14213E-03

rezultat of identification ND= 3

A1= -.35703 A2= .41887E-01 A3= -.10388E-02

B1= -.17768E-03 B2= 4.9999 B3= -.48533

C0= .19696E-03

CKO= .13943E-03

rezultat of identification ND= 4

A1= -.10243E-01 A2= -.61121E-01 A3= .80717E-02

B1= 4.9998 B2= 1.2485 B3= -.64322E-01

C0= .30525E-03

CKO= .13293E-03

rezultat of identification ND= 5

A1= 646.76 A2= -172.26 A3= 13.158

B1= 3233.9 B2= -21.478 B3= 6.1346

C0= -.19609

CKO= .74519

rezultat of identification ND= 6

A1= -.56044E-01 A2= 330.98 A3= -54.860

B1= 1653.6 B2= 155.01 B3= 12.592

C0= -.79423E-01

CKO= .76317

rezultat of identification ND= 7

A1= .14801E-01 A2= .20602 A3= -23.377

B1= -117.88 B2= -30.870 B3= -6.2714

C0= -.21560E-02

CKO= .78029

rezultat of identification ND= 8

A1= .12913 A2= .34425 A3= .18700

B1= -1.2883 B2= -1.6156 B3= .32495

C0= .41111

CKO= .69853

CKOmin=CKO(4)= .13293E-03

Optimal model: ND= 4

A1= -.10243E-01 A2= -.61121E-01 A3= .80717E-02

B1= 4.9998 B2= 1.2485 B3= -.64322E-01

C0= .30525E-03

CKO= .13293E-03

ishodn. data for identification

L2=31

Output for model optimal

INPUT OUTPUT

-1.4710 5.1980

-.5785 6.3200

1.0460 5.8260

-.5766 3.3380

-.4351 -3.6970

-1.2630 -8.3750

.7064 -5.0110

.9887 4.0708

.0995 -1.7370

.8031 -2.6910

.3668 -6.9860

-1.1940 1.7610

.1346 5.5190

.6195 1.9070

.0273 4.4190

-.1278 2.9480

-.3297 -5.2780

-.8129 -.7508

.6869 2.9880

1.1880 .9287

-1.2950 -.4464

.5024 -1.7820

-.6517 -4.5210

-.2448 2.2890

-.7592 6.6120

.5688 -4.7900

-.9917 1.1560

-1.4290 -2.8820

.2751 -1.9900

-.6137 -4.2660

-.3316 1.7700

-.3619 -4.4256

-.1892 -8.3219

-.0478 -.7146

.3277 -3.1129

.9536 -2.4499

.8451 -2.3932

.2682 -1.5253

.9878 -.5938

1.5210 1.5112

.3804 5.1718

.3346 5.5453

-1.3520 2.6957

.3365 5.5444

-1.5840 8.9978

.1512 4.1470

.3082 2.5980

.4029 -6.1587

-.0488 .0355

-.2819 -7.8093

-2.0330 -1.2711

.7191 1.3413

.9858 2.3889

.0341 .3562

.6283 -1.3571

-.7133 -10.5245

.5209 .8819

.1005 5.3344

-.7103 1.5488

.1024 3.4556

Stop - Program terminated.

Квазирекуррентный алгоритм идентификации .

Рекуррентные алгоритмы наименьших квадратов характеризуются малым объемом вычислений благодаря использованию вектора коррекции γ _i и обратной матрицы (Ф^ТФ), однако по этой же причине излишняя сложность формируемых моделей не влияет на работоспособность рекуррентных алгоритмов и не может быть обнаружена.

Достоинствами нерекуррентного и рекуррентного алгоритмов идентификации обладает квазирекуррентный алгоритм идентификации.

В основе квазирекуррентного алгоритма идентификации лежит то, что в уравнении =[Ф^ТФ]^-1Ф^ТХ, описывающем нерекуррентный алгоритм идентификации, сомножители Ф^ТФ и Ф^ТХ имеют постоянные размерности 2 m x 2 m и m и могут быть определены рекуррентно.

[Ф^ТФ] _k - [Ф^ТФ] _k-1 + φ( k) φ^Т( k), [Ф^ТФ]₀ = 0;

[Ф^ТХ] _k - [Ф^ТХ] _k-1 + x( k) φ( k), [Ф^ТХ]₀ = 0.

Используя приведенные выше соотношения, после введения новой матрицы D( k) и вектора C( k), определяемых уравнениями:

D( k)=[Ф^ТФ] _k;

C( k)=[Ф^ТХ] _k,

можно записать соотношения, определяющие квазирекуррентный алгоритм идентификации:

D(k)= D(k-1)+ φ(k)φ^Т (k), D(0)=0;

C( k)= C( k-1)+ x( k) φ( k), C(0)=0,

где C( k) - вектор размерности 2 m; D( k) - матрица размерности 2 m x 2 m; x( k) - выходная координата; и уравнение для оценки вектора параметров:

(3.1)

После вычисления оценки вектора параметров на каждом шаге можно оценить качество используемой модели динамического объекта по значению функции потерь:

(3.2)

где S( k) - скаляр, рекуррентно вычисляемый по формуле

S( k) = λS( k-1) + x²( k), S(0) = 0.

Матрица D( k) является симметричной и при программной реализации алгоритма во всех вычислениях может использоваться либо ее верхняя, либо нижняя треугольная часть. Затраты времени и памяти ЦВМ на реализацию квазирекуррентного алгоритма наименьших квадратов близки к затратам на реализацию рекуррентных алгоритмов, а сами алгоритмы совместимы в одной программе, однако при этом квазирекуррентный алгоритм позволяет оценивать степень соответствия сложности модели имеющимся экспериментальным данным.

Двух ступенчатость квазирекуррентного алгоритма наименьших квадратов позволяет производить по результатам обработки экспериментальных данных с помощью уравнений (3.1), (3.2) поиск линейных моделей оптимальной сложности на основании сформированных матриц C( k) и D( k).

Эта возможность обусловлена тем, что матрицы C( k) и D( k) состоят из клеточных матриц, в каждой из которых хранится информация о составляющих вектора параметров. Пусть общая модель имеет пять входов и один выход, тогда структура матриц C( k) и D( k) имеет вид, представленный на рис. 3.7.

Для поиска моделей оптимальной сложности формируются матрицы С' и D'. В качестве примера на рис. 3.8 показаны матрицы С' и D', сформированные для идентификации одной из частных моделей вида „2 входа -1 выход".

Размерность соответствующих клеточных матриц, используемых для формирования матриц D' и С', может быть равной или меньшей, чем в матрицах D_k и C_k. Вырожденность матрицы D' свидетельствует о чрезмерной сложности анализируемой модели; если матрица D' невырождена, то после вычисления соответствующей частной модели оценки вектора параметров по формуле:

можно оценить качество сформированной модели по соответствующему ей значению функции потерь:

D_k C_k

Й вход 1-й вход

2-й вход 1 1^’ 2^’ 4^’ 2-й вход 7

Й вход 3-й вход

4-й вход 2 3 3^’ 5^’ 4-й вход 8

Й вход 5-й вход

Выход 4 5 6 6^’ Выход 9

1-й вход

2-й вход

3-й вход

4-й вход

5-й вход

Выход

Рис. 3.7.

1^’ 1	2^’	4^’
2^’	1^’ 1	5^’
4^’	5^’	1^’ 1

D ^’ C ^’

1-й вход 7

4-й вход 8

Выход 9

2-й вход

4-й вход

Выход

Рис.3.8.

Поскольку число частных моделей у сложного динамического объекта может быть значительным, для поиска модели объекта оптимальной сложности целесообразно использовать специальный алгоритм селекции моделей, например метод группового учета аргументов (МГУА).

Квазирекуррентный алгоритм идентификации позволяет вычислять вектор как в темпе со временем на каждом шаге обработки данных, так и в произвольные моменты времени. При использовании квазирекуррентного алгоритма в темпе со временем вычислительные затраты примерно в 1,5 раза больше, чем при использовании рекуррентных алгоритмов. Однако во многих задачах не требуется вычислять оценку вектора с той же частотой, с какой поступает информация. А уже при частоте оценки втрое меньшей частоты поступления данных вычислительные затраты квазирекуррентного алгоритма и рекуррентных алгоритмов примерно равны.

Объем памяти ЭВМ, необходимый для работы квазирекуррентного алгоритма при размерности вектора параметров , равной 41, примерно в 2,5 раза больше, чем для рекуррентного алгоритма, и примерно в 10 раз меньше, чем для нерекуррентного алгоритма.

Дата добавления: 2018-11-24; просмотров: 73; Мы поможем в написании вашей работы!

Поделиться с друзьями:

Мы поможем в написании ваших работ!