SUBROUTINE AOLSYS(SYSFIL,ANLFIL,ERRCD,ERRMSG) C C FUNCTION: CF CF C USAGE: CU CU C INPUTS: CI CI C OUTPUTS: CO CO ERROR_CODE integer CO value range description CO INSYS code [0..100] see INSYS CO 100 + rg code [101..200] see rg (eispack) CO 200 - mvzero code [201..300] see mvzero CO 300 + ahcon code [301..400] see ahcon (CC calc) CO 400 + ahcon code [401..500] see ahcon (CO calc) CO C ALGORITHM: CA CA C MACHINE DEPENDENCIES: CM CM C HISTORY: CH CH written by: J. Douglas Birdwell CH date: 23-july-1986 CH current version: 2.1 CH modifications: standardized,modified to accomodate CH header.c bb: 8-86. CH fixed problem with z2 and z1 CH bb and jdb :8-86. CH added dpcom: 7/16/88 jdb CH C ROUTINES CALLED: CC CC C COMMON MEMORY USED: CM CM DPCOM -- see dpcommon.f and dpcom.f CM C---------------------------------------------------------------------- C written for: The CASCADE Project C Oak Ridge National Laboratory C U.S. Department of Energy C contract number DE-AC05-840R21400 C subcontract number ########## C organization: xxxxxx C---------------------------------------------------------------------- C THIS SOFTWARE IS IN THE PUBLIC DOMAIN C NO RESTRICTIONS ON ITS USE ARE IMPLIED C---------------------------------------------------------------------- C INCLUDE 'Parameter.f' C DOUBLE PRECISION A(SIZE,SIZE) DOUBLE PRECISION B(SIZE,SIZE) DOUBLE PRECISION C(SIZE,SIZE) DOUBLE PRECISION D(SIZE,SIZE) DOUBLE PRECISION ATEMP(SIZE,SIZE) DOUBLE PRECISION Z1(SIZE2,SIZE2) DOUBLE PRECISION Z2(SIZE2,SIZE2) DOUBLE PRECISION REZERO(SIZE) DOUBLE PRECISION IMZERO(SIZE) DOUBLE PRECISION REPOLE(SIZE) DOUBLE PRECISION IMPOLE(SIZE) DOUBLE PRECISION POLES(SIZE,2) DOUBLE PRECISION ZEROS(SIZE,2) DOUBLE PRECISION BETA(SIZE) DOUBLE PRECISION SUM(SIZE) DOUBLE PRECISION DUMMY(SIZE) DOUBLE PRECISION SIGMA(SIZE) DOUBLE PRECISION U DOUBLE PRECISION V DOUBLE PRECISION EIGVEC(SIZE,SIZE) DOUBLE PRECISION EPS DOUBLE PRECISION EPS1 DOUBLE PRECISION EPSLON DOUBLE PRECISION WNZERO DOUBLE PRECISION THZERO DOUBLE PRECISION FREQ DOUBLE PRECISION WNMIN DOUBLE PRECISION WNMAX C INTEGER NZEROS INTEGER NRANK INTEGER NSTATS INTEGER NINPS INTEGER NOUTS INTEGER ERRCD INTEGER IPVT(SIZE) INTEGER JPVT(SIZE) INTEGER NPORHP INTEGER NPJWAX INTEGER NPORIG INTEGER NPOLHP INTEGER NZORHP INTEGER NZJWAX INTEGER NZORIG INTEGER NZOLHP INTEGER NOBSMD INTEGER NCONMD INTEGER IDXSIG(SIZE) INTEGER NCON0(SIZE) INTEGER NOBS0(SIZE) INTEGER IDXPOR(SIZE) INTEGER ISEED INTEGER LENGTH INTEGER MERR INTEGER MAXLMN C LOGICAL INSBSP LOGICAL CON(SIZE) C CHARACTER*(*) ANLFIL CHARACTER*(*) SYSFIL CHARACTER*(*) ERRMSG CHARACTER*30 STRING C INCLUDE 'dpcom.f' C ERRCD=0 C C READ IN THE SYSTEM FROM THE SYSTEM CONTAINER FILE C CALL INSYS(SYSFIL,NINPS,NOUTS,NSTATS, 2 SIZE,SIZE,SIZE,A,B,C,D,ERRCD) CLOSE (UNIT=UNIT1) IF (ERRCD .NE. 0) THEN ERRMSG = 'AOLSYS: Fatal error from INSYS '// 1 'accessing '//SYSFIL RETURN END IF C C--calculate the poles of the system C CALL SAVE (SIZE,SIZE,NSTATS,NSTATS,A,ATEMP) CALL RG (SIZE,NSTATS,ATEMP,REPOLE,IMPOLE, . 0,U,IPVT,DUMMY,ERRCD) C IF (ERRCD .NE. 0) THEN ERRMSG='AOLSYS: Fatal error from '// 2 'eispack (RG) ERRCD = 100 + return code: ' ERRCD = 100 + ERRCD RETURN END IF C C--calculate multivariable zeros of system C CALL MVZERO (NOUTS,NINPS,NSTATS, 2 SIZE,SIZE,SIZE2, 3 A,B,C,D,Z1,Z2,NZEROS,NRANK, 4 REZERO,IMZERO, 5 BETA,SUM,DUMMY) C IF (NZEROS .LT. 0) THEN ERRMSG = ' AOLSYS: Zero calculation failed.' ERRCD = 200 - NZEROS RETURN END IF C C--find maximum natural frequency C WNMAX = 0.D0 C DO 10, I = 1, NSTATS C FREQ = SQRT(REPOLE(I)**2 2 + IMPOLE(I)**2) IF (FREQ .GT. WNMAX) WNMAX = FREQ C 10 CONTINUE C C--compute machine precision C EPS = 1.D0 20 CONTINUE EPS = EPS / 2.D0 EPS1 = EPS + 1.D0 IF (EPS1 .GT. 1.D0) GO TO 20 EPS = 2.D0 * EPS C C--set width of "jw" axis C EPSLON = EPS**0.3333333333333333D0 * MAX(WNMAX,1.D0) C C--find minimum natural frequency C WNMIN = 1.D30 C DO 30, I = 1, NSTATS C FREQ = SQRT(REPOLE(I)**2 2 + IMPOLE(I)**2) IF (FREQ .LT. WNMIN .AND. FREQ .GT. EPSLON) WNMIN = FREQ C 30 CONTINUE C IF (WNMIN .EQ. 1.D30) WNMIN = 0.001D0 C C--count zeros in ORHP, on jw axis, at 0, and in OLHP. C NZORHP = 0 NZJWAX = 0 NZORIG = 0 NZOLHP = 0 C DO 40, I = 1, NZEROS C C Save the zeros of the system C ZEROS(I,1) = REZERO(I) ZEROS(I,2) = IMZERO(I) C IF (REZERO(I) .LT. -EPSLON) THEN NZOLHP = NZOLHP + 1 ELSE IF (REZERO(I) .GT. EPSLON) THEN NZORHP = NZORHP + 1 ELSE IF (DABS(IMZERO(I)) .GT. EPSLON) THEN NZJWAX = NZJWAX + 1 ELSE NZORIG = NZORIG + 1 END IF C 40 CONTINUE C C--count poles in ORHP, on jw axis, at 0, and in OLHP. C NPORHP = 0 NPJWAX = 0 NPORIG = 0 NPOLHP = 0 C DO 50, I = 1, NSTATS C C Save the poles of the system. C POLES(I,1) = REPOLE(I) POLES(I,2) = IMPOLE(I) C IF (REPOLE(I) .LT. -EPSLON) THEN NPOLHP = NPOLHP + 1 ELSE IF (REPOLE(I) .GT. EPSLON) THEN NPORHP = NPORHP + 1 ELSE IF (DABS(IMPOLE(I)) .GT. EPSLON) THEN NPJWAX = NPJWAX + 1 ELSE NPORIG = NPORIG + 1 C C--keep track of where the eigenvectors corresponding to zero eigenvalues C--are for controllability calculations C IDXPOR(NPORIG) = I C END IF C 50 CONTINUE C C--create output file C OPEN (UNIT=UNIT1,FILE=ANLFIL) C C--write pole data C WRITE (UNIT1,*) NPORHP, 2 NPJWAX, 3 NPORIG, 4 NPOLHP C C--write zero data C WRITE (UNIT1,*) NZORHP, 2 NZJWAX, 3 NZORIG, 4 NZOLHP C C--find natural frequency of smallest non-minimum phase zero C IF (NZORHP .EQ. 0) THEN WNZERO = 0.D0 ELSE WNZERO = 1.0D30 DO 60, I = 1, NZEROS IF (REZERO(I) .GT. EPSLON) THEN THZERO = SQRT(REZERO(I)**2 2 + IMZERO(I)**2) WNZERO = MIN(WNZERO,THZERO) END IF 60 CONTINUE END IF C WRITE (UNIT1,*) WNZERO C C--begin controllability and observability calculations C ISEED = 0 CALL AHCON (SIZE,NSTATS,NINPS,A,B,REPOLE, . IMPOLE,REZERO,IMZERO, . DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD) IF (ERRCD .NE. 0) THEN ERRCD = 300 + ERRCD ERRMSG = ' AOLSYS: Controllability calculation failed.'// 1 ' System is not controllable.' RETURN END IF NCONMD = 0 DO 70, I = 1, NSTATS IF (CON(I)) NCONMD = NCONMD + 1 70 CONTINUE C CALL TRNATB (SIZE,SIZE2,NSTATS,NSTATS,A,Z1) CALL TRNATB (SIZE,SIZE2,NOUTS,NSTATS,C,Z2) CALL AHCON (SIZE2,NSTATS,NOUTS,Z1,Z2,REPOLE, . IMPOLE,REZERO,IMZERO, . DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD) IF (ERRCD .NE. 0) THEN ERRCD = 400 + ERRCD ERRMSG = ' AOLSYS: Observability calculation failed.'// 1 ' System is not observable.' RETURN END IF NOBSMD = 0 DO 80, I = 1, NSTATS IF (CON(I)) NOBSMD = NOBSMD + 1 80 CONTINUE C C--now compute the number of controllable s=0 modes from each input C C--if there are no poles at s=0, skip the work C IF (NPORIG .GT. 0) THEN C C--for each input, check controllability using ahcon C DO 90, I = 1, NINPS C CALL AHCON (SIZE,NSTATS,1,A,B(1,I),REPOLE, . IMPOLE,REZERO,IMZERO, . DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD) IF (ERRCD .NE. 0) THEN CALL ICNVRT(0,I,STRING,LENGTH,MERR) ERRCD = 500 + ERRCD ERRMSG=' AOLSYS: Controllability calculation '// . 'failed for input '//STRING(1:LENGTH)//'.' RETURN END IF C C--check for each zero eigenvalue C NCON0(I) = 0 DO 100, J = 1, NSTATS FREQ = SQRT(REPOLE(J)**2 . + IMPOLE(J)**2) C IF (FREQ .LT. EPSLON) THEN C C--this is a zero eigenvalue C--now, look for the eigenvalue closest to this eigenvalue C DO 110, K = 1, NSTATS IF (IPVT(K) .EQ. J) INDEX = K 110 CONTINUE C C--and check to see if it is sufficiently far away from the 0 eigenvalue C IF (CON(INDEX)) NCON0(I) = NCON0(I) + 1 END IF 100 CONTINUE 90 CONTINUE C C--else, there are no s=0 modes and all this can be skipped C ELSE C DO 120, I = 1, NINPS NCON0(I) = 0 120 CONTINUE C END IF C C--now compute the number of observable s=0 modes from each output C C--if there are no poles at s=0, skip the work C IF (NPORIG .GT. 0) THEN C C--remember that z1 contains a-transpose and z2 contains b-transpose C C--for each output, check observability using ahcon C DO 130 I = 1, NOUTS CALL AHCON (SIZE2,NSTATS,1,Z1,Z2(1,I),REPOLE, . IMPOLE,REZERO,IMZERO, . DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD) IF (ERRCD .NE. 0) THEN CALL ICNVRT(0,I,STRING,LENGTH,MERR) ERRCD = 600 + ERRCD ERRMSG = ' AOLSYS: observability calculation '// . 'failed for output '//STRING(1:LENGTH)//'.' RETURN END IF C C--check for each zero eigenvalue C NOBS0(I) = 0 DO 140, J = 1, NSTATS FREQ = SQRT(REPOLE(J)**2 . + IMPOLE(J)**2) C IF (FREQ .LT. EPSLON) THEN C C--this is a zero eigenvalue C--now, look for the eigenvalue closest to this eigenvalue C DO 150, K = 1, NSTATS IF (IPVT(K) .EQ. J) INDEX = K 150 CONTINUE C C--and check to see if it is sufficiently far away from the 0 eigenvalue C IF (CON(INDEX)) NOBS0(I) = NOBS0(I) + 1 END IF 140 CONTINUE 130 CONTINUE C C--else, there are no s=0 modes and all this can be skipped C ELSE C DO 160, I = 1, NOUTS NOBS0(I) = 0 160 CONTINUE C END IF C C--print out number of controllable modes of the system C WRITE (UNIT1,*) NCONMD C C--print out number of observable modes of the system C WRITE (UNIT1,*) NOBSMD C C--PRINT OUT VECTOR OF NUMBER OF CONTROLLABLE S=0 MODES OF C--the system from input i, for i=1,..,NINPS C DO 170, I = 1, NINPS WRITE (UNIT1,*) NCON0(I) 170 CONTINUE C C--print out vector of number of observable s=0 modes of C--the system from output i, for i =1,..,NOUTS C DO 180, I = 1, NOUTS WRITE (UNIT1,*) NOBS0(I) 180 CONTINUE C WRITE (UNIT1,*) WNMIN, WNMAX C C Write out poles of system C WRITE(UNIT1,185) 185 FORMAT(1X,'POLES:') C DO 190, I=1,NSTATS WRITE(UNIT1,*) POLES(I,1),POLES(I,2) 190 CONTINUE C C Write out zeros of system C WRITE(UNIT1,195) 195 FORMAT(1X,'ZEROS:') C DO 200, I=1,NZEROS WRITE(UNIT1,*) ZEROS(I,1),ZEROS(I,2) 200 CONTINUE C C--end of processing C CLOSE (UNIT=UNIT1) RETURN END