      SUBROUTINE INITAL
C...
C...  PI CONTROL OF A BINARY BATCH DISTILLATION COLUMN
C...
C...  THE FOLLOWING EQUATIONS MODEL A BATCH DISTILLATION COLUMN WITH
C...  PI CONTROL OF THE DISTILLATE RATE
C...
C...  KETTLE       HB*DXB/DT = R*X1 - V*YB, XB(0) = XB0
C...
C...  PLATE 1      H*DX1/DT = R*X2 + V*YB - R*X1 - V*Y1, X1(0) = X10
C...
C...  PLATE 2      H*DX2/DT = R*X3 + V*Y1 - R*X2 - V*Y2, X2(0) = X20
C...
C...  PLATE N      H*DXN/DT = R*XD + V*YN-1 - R*XN - V*YN, XN(0) = XN0
C...
C...  CONDENSER    HC*DXD/DT = V*YN - R*XD - D*XD, XD(0) = XD0
C...
C...  EQUILIBRIUM  YN = A*XN/(1 + (A - 1)*XN)
C...
C...                                        T
C...  CONTROLLER   Z = ZS + KC*(E + (1/TI)*INT E*DT)
C...                                        0
C...               E = XD - XDSET
C...
C...               D = CV*Z*(HC**0.5)
C...
C...               R = V - D
C...
C...  WHERE
C...
C...     X1,X2,    LIQUID COMPOSITIONS ON TRAYS 1, 2 ..., N
C...     ...,XN
C...
C...     Y1,Y2,    VAPOR COMPOSITIONS ABOVE TRAYS 1, 2,..., N
C...     ...,YN
C...
C...     A         RELATIVE VOLATILITY
C...
C...     XB,XD     REBOILER, CONDENSER LIQUID COMPOSITIONS, RESPECTIVELY
C...
C...     XDSET     CONDENSER LIQUID COMPOSITION SET POINT
C...
C...     YB        COMPOSITION OF VAPOR FROM THE REBOILER
C...
C...     T         TIME
C...
C...     HB,H,HC   REBOILER HOLDUP, HOLDUPS ON TRAYS 1 TO N, CONDENSER
C...               HOLDUP, RESPECTIVELY
C...
C...     V         VAPOR BOILUP RATE
C...
C...     R         REFLUX RATE (LIQUID RATE)
C...
C...     D         DISTILLATE RATE
C...
C...     CV        CONTROL VALVE CONSTANT
C...
C...     Z         VALVE STEM POSITION
C...
C...     ZS        STEADY STATE VALVE STEM POSITION
C...
C...     KC        CONTROLLER GAIN
C...
C...     TI        CONTROLLER INTEGRAL TIME
C...
C...     E         ERROR IN DISTILLATE COMPOSITION = XD - XDSET
C...
C...  THE PROBLEM IS BASICALLY TO SELECT THE CONTROLLER GAIN AND
C...  INTEGRAL TIME TO MAINTAIN SATISFACTORY CONTROL OF THE DISTILLATE
C...  COMPOSITION.  TWO FIGURES OF MERIT FOR THE CONTROLLER PERFORMANCE
C...  ARE TO BE COMPUTED
C...
C...           T              T
C...     I1 = INT D*DT, I2 = INT (XD*D)*DT
C...           0              0
C...
C...  INTEGRAL I2 GIVES THE TOTAL YIELD OF THE LIGHT COMPONENT (WITH
C...  MOL FRACTION XD IN THE DISTLLATE), AND THE RATIO I2/I1 GIVES THE
C...  TIME AVERAGE DISTILLATE COMPOSITION.  THESE INTEGRALS CAN BE
C...  EVALUATED BY PROGRAMMING THEM AS DIFFERENTIAL EQUATIONS, I.E.,
C...
C...    DI1/DT = D, I1(0) = 0, DI2/DT = XD*D, I2(0) = 0
C...
      COMMON/T/  TIME, NSTOP, NORUN
     1      /Y/    XD, X(10),    XB,    IE,    I1,    I2
     2      /F/   XDT,XT(10),   XBT, DIEDT, DI1DT, DI2DT
     3      /V/     A,     V,     R,     H,    HC,    HB,    YB,
     4          Y(10),    CV,    ZS,     Z,     E,    KC,    TI,
     5          XDSET, RATIO,     D,     N
      REAL         IE,    I1,    I2,    KC
C...
C...  SET THE MODEL PARAMETERS
      N=10
      A=2.5
      V=5.
      H=0.5
      HC=3.
      HB=100.
      CV=2.309
      ZS=0.5
      XDSET=0.85
C...
C...  SET THE CONTROLLER CONSTANTS
      IF(NORUN.EQ.1)THEN
      KC=0.
      TI=1.0E+10
      ELSE IF(NORUN.EQ.2)THEN
      KC=10.
      TI=25.
      END IF
C...
C...  SET THE INITIAL CONDITIONS
10    XD=0.5
      DO 11 I=1,N
      X(I)=0.5
11    CONTINUE
      XB=0.5
      IE=0.
      I1=0.
      I2=0.
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/  TIME, NSTOP, NORUN
     1      /Y/    XD, X(10),    XB,    IE,    I1,    I2
     2      /F/   XDT,XT(10),   XBT, DIEDT, DI1DT, DI2DT
     3      /V/     A,     V,     R,     H,    HC,    HB,    YB,
     4          Y(10),    CV,    ZS,     Z,     E,    KC,    TI,
     5          XDSET, RATIO,     D,     N
      REAL         IE,    I1,    I2,    KC
C...
C...  MODEL ALGEBRA
C...
C...  VAPOR COMPOSITIONS
      DO 1 I=1,N
      Y(I)=EQUIL(A,X(I))
1     CONTINUE
C...
C...  CONTROLLER EQUATIONS
      YB=EQUIL(A,XB)
      E=XD-XDSET
      Z=ZS+KC*(E+(1./TI)*IE)
      IF(Z.LT.0.)Z=0.
      IF(Z.GT.1.)Z=1.
      D=CV*Z*SQRT(HC)
      R=V-D
C...
C...  CONDENSER
      XDT=(1./HC)*(V*Y(N)-R*XD-D*XD)
C...
C...  PLATE N
      XT(N)=(1./H)*(R*XD+V*Y(N-1)-R*X(N)-V*Y(N))
C...
C...  PLATES 2 TO N-1
      NM1=N-1
      DO 2 I=2,NM1
      XT(I)=(1./H)*(R*X(I+1)+V*Y(I-1)-R*X(I)-V*Y(I))
2     CONTINUE
C...
C...  PLATE 1
      XT(1)=(1./H)*(R*X(2)+V*YB-R*X(1)-V*Y(1))
C...
C...  KETTLE
      XBT=(1./HB)*(R*X(1)-V*YB)
C...
C...  TIME INTEGRALS
      DIEDT=E
      DI1DT=D
      DI2DT=XD*D
      RETURN
      END
      FUNCTION EQUIL(A,X)
C...
C...  FUNCTION EQUIL CALCULATES THE EQUILIBRIUM VAPOR COMPOSITION
C...  FROM ROAULT*S LAW
C...
C...  ARGUMENT LIST
C...
C...     A     RELATIVE VOLATILITY (INPUT)
C...
C...     X     LIQUID MOL FRACTION, FOR WHICH THE EQUILIBRIUM VAPOR
C...           COMPOSITION IS TO BE COMPUTED (INPUT)
C...
      EQUIL=A*X/(1.+(A-1.)*X)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/  TIME, NSTOP, NORUN
     1      /Y/    XD, X(10),    XB,    IE,    I1,    I2
     2      /F/   XDT,XT(10),   XBT, DIEDT, DI1DT, DI2DT
     3      /V/     A,     V,     R,     H,    HC,    HB,    YB,
     4          Y(10),    CV,    ZS,     Z,     E,    KC,    TI,
     5          XDSET, RATIO,     D,     N
      REAL         IE,    I1,    I2,    KC
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION PT(51), PXD(51), PD(51), PZ(51), PE(51), PI1(51),
     1         PI2(51),  PR(51)
      DATA IP/0/
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)WRITE(NO,1)
1     FORMAT(9X,1HT,10X,2HXD,11X,1HD,11X,1HZ,11X,1HE,10X,2HIE,
     1      10X,2HI1,10X,2HI2,7X,5HI2/I1)
C...
C...  RPINT THE NUMERICAL SOLUTION.  NOTE - THE RATIO I2/I1 CAN BE
C...  INDETERMINATE (0/0).  IF THIS OCCURS, THIS RATIO IS ARBITRARILY
C...  SET TO ZERO
      IF((I1.EQ.0.).AND.(I2.EQ.0.))RATIO=0.
      IF((I1.EQ.0.).AND.(I2.NE.0.))STOP
      IF((I1.NE.0.).AND.(I2.NE.0.))RATIO=I2/I1
      WRITE(NO,2)TIME,XD,D,Z,E,IE,I1,I2,RATIO
2     FORMAT(F10.2,8F12.3)
C...
C...  STORE THE SOLUTION FOR PLOTTING
       IP=IP+1
       PT(IP)=TIME
      PXD(IP)=XD
       PD(IP)=D
       PZ(IP)=Z
       PE(IP)=E
      PI1(IP)=I1
      PI2(IP)=I2
       PR(IP)=RATIO
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.51)RETURN
      CALL SPLOTS(1,IP,PT,PXD)
      WRITE(NO,10)
10    FORMAT(1H ,//,8H XD VS T)
      CALL SPLOTS(1,IP,PT,PD )
      WRITE(NO,11)
11    FORMAT(1H ,//,8H  D VS T)
      CALL SPLOTS(1,IP,PT,PZ )
      WRITE(NO,12)
12    FORMAT(1H ,//,8H  Z VS T)
      CALL SPLOTS(1,IP,PT,PE )
      WRITE(NO,13)
13    FORMAT(1H ,//,8H  E VS T)
      CALL SPLOTS(1,IP,PT,PI1)
      WRITE(NO,14)
14    FORMAT(1H ,//,8H I1 VS T)
      CALL SPLOTS(1,IP,PT,PI2)
      WRITE(NO,15)
15    FORMAT(1H ,//,8H I2 VS T)
      CALL SPLOTS(1,IP,PT,PR )
      WRITE(NO,16)
16    FORMAT(1H ,//,11H I2/I1 VS T)
      IP=0
      RETURN
      END
