optexpsun.f

对工业生产过程结晶过程的一个仿真程序软件包

* Copyright c 1998-2002 The Board of Trustees of the University of Illinois
* 		  All rights reserved.
* Developed by:	Large Scale Systems Research Laboratory
*               Professor Richard Braatz, Director
*               Department of Chemical Engineering
*		University of Illinois
*		http://brahms.scs.uiuc.edu
* 
* Permission hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal with the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the 
* Software is furnished to do so, subject to the following conditions:
* 		1. Redistributions of source code must retain the above copyright
*		   notice, this list of conditions and the following disclaimers.
*		2. Redistributions in binary form must reproduce the above 
*		   copyright notice, this list of conditions and the following 
*		   disclaimers in the documentation and/or other materials 
*		   provided with the distribution.
*		3. Neither the names of Large Scale Research Systems Laboratory,
*		   University of Illinois, nor the names of its contributors may
*		   be used to endorse or promote products derived from this 
*		   Software without specific prior written permission.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
* THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
* OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*
*	optexpsun.f
*	
*
*	This program calculates the optimum temperature profile
*	and initial seed distribution for optimum experimental 
*	design.  The temperature profile is formed by piecewise 
*	linear trajectories (See the subroutine Temp).  The 
*	following parameters (in the given subroutine) should be 
*	set to the number of discretizations (i.e. the  number of 
*	linear trajectories):
*
*		Subroutine/Program	Parameter
*		------------------	---------
*		Main			ntemp
*		FCN			Ntemp1
*		cntr			Ntemp2
*		Temp			Ntemp3 
*
*	The the initial seed distribution is characterized
*	by the total mass, mean size, and width of the distribution.
*
*	Parameter inputs are growth and nucleation 
*       kinetic parameters (g, kg, b, and kb)
*
*	The optimization problem is solved using the sequential
*       quadratic program subroutine FFSQP by Jian L. Zhou, Andre L. 
*	Tits, and C.T. Lawrence.  FFSQP and the attached subroutines
*	are given below.
*	
*
*       Date:    July 6, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign
*
************************************************************************
      PROGRAM MAIN

*	The main program initializes the variables used by the
*	FFSQP subroutine.  Explanation of the variables is given 
*	in the FFSQP subroutine, except of ntemp, which is the number of
*	temperature discretization.  After initialization the program
*	calls the FFSQP subroutine to solve the nonlinear constrained 
*	optimization problem.

      INTEGER nparam, nf, nineqn, nineq, neqn, neq, iwsize, nwsize
      INTEGER mode, iprint,miter
      INTEGER inform, ntemp
      PARAMETER(ntemp = 8, nparam=ntemp+3,nf=1, 
     &		nineq=1, neq=0)
      PARAMETER(mode=100,miter=10000)
      PARAMETER(iwsize=6*nparam+8*(nineq+neq)+7*(nf)+30) 
      PARAMETER(nwsize=4*nparam**2+5*(nineq+neq)*nparam+
     &		3*(nf)*nparam+
     &		26*(nparam+nf)+45*(nineq+neq)+100) 
      INTEGER iw(iwsize)
      INTEGER I
      REAL*8 bigbnd, eps, epseqn, udelta
      REAL*8 x(nparam), bl(nparam), bu(nparam)
      REAL*8 f(nf),g(nineq+neq),w(nwsize)
      REAL*8 scale(nparam)
      EXTERNAL FCN,cntr,grobfd,grcnfd

      bigbnd=1.0D12
      eps=1.0D-8
      epseqn=0.0D0
      udelta=0.0D0
      iprint=3
      nineqn=0
      neqn=0

* Initial guess:
* x(1) to x(ntemp) are the slopes of the linear pieces for the
* the temperature profile.  x(ntemp+1) is the total mass of seed
* in grams. x(ntemp+2) is the first moment of the seed distribution 
* in micron/g solvent.  x(ntemp+3) is the width (in microns) of the 
* domain of the crystal size distribution function.

      x(1) = -1.0000000000000D-01
      x(2) =  2.6636175466309D-33
      x(3) = -1.0000000000000D-01
      x(4) = -1.0000000000000D-01
      x(5) = -1.0000000000000D-01
      x(6) = -1.0000000000000D-01
      x(7) = 4.0325592487296D-34
      x(8) = 4.0325592487296D-34
      x(9) = 5.0D0
      x(10) = 600.0D0
      x(11) = 95.0D0

*Scaling factors for the parameters

      DO 202 I = 1, ntemp
202	scale(I)=1.0D0
        scale(ntemp+1)=1.0D-3
        scale(ntemp+2)=1.0D-3
        scale(ntemp+3)=1.0D0

      DO 201 I = 1, nparam
201	x(I)=scale(I)*x(I)


* Lower bound for the parameter:
      DO 666 I = 1, ntemp
666	bl(I)=-0.1D0
      bl(ntemp+1)=5.0D0*scale(ntemp+1)
      bl(ntemp+2)=5.0D0*scale(ntemp+2)
      bl(ntemp+3)=5.0D0*scale(ntemp+3)

* Upper bound for the parameters:
      DO 667 I = 1, ntemp
667	bu(I)=0.0D0
      bu(ntemp+1)=110000.0D0*scale(ntemp+1)
      bu(ntemp+2)=600.0D0*scale(ntemp+2)
      bu(ntemp+3)=95.0D0*scale(ntemp+3)

      call FFSQP(nparam,nf,nineqn,nineq,neqn,neq,mode,iprint,
     *           miter,inform,bigbnd,eps,epseqn,udelta,bl,bu,x,f,g,
     *           iw,iwsize,w,nwsize,FCN,cntr,grobfd,grcnfd)

      PRINT*,'Final Solution'
      DO 301 I = 1, nparam 
301      	PRINT*,x(I)/scale(I)
      PRINT*,'objective = ', f


      STOP
      END

******************************************************************
*
*       SUBROUTINE FCN(N,j,x,fj)
*
*       N = number of parameters to be determined
*       x = vector of length N, parameters to be determined
*
*
*       This subroutine simulates the operation of an industrial
*       crystallizer scaled-up from the experimental batch cooling
*       crystallizer described in S. M. Miller's Ph.D. thesis
*       published at the University of Texas at Austin in 1993.
*
*
*	The following subroutines are required for FCN:
*	
*		Temp
*		MOMENTS
*		MOMENTSJ
*		Csat
*		growth
*		birth
*

      SUBROUTINE FCN(Nvar,jjj,x,fj)

      INTEGER Nvar,jjj,Ntemp1
      INTEGER NSTEP, NEQ, MXPARM, NTHETA, NU
      PARAMETER (NSTEP=161,NEQ=33,MXPARM=50,NTHETA=4, NU=6)
      PARAMETER (Ntemp1 = 8)
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
      INTEGER I, J, K, M, NNN
      REAL*8 x(*), fj
      REAL*8 Coeff(Ntemp1+3)
      REAL*8 T, Y(NEQ)
      REAL*8 F(2*(NSTEP-1), NTHETA), FTVF(NTHETA, NTHETA)
      REAL*8 FTV(NTHETA,2*(NSTEP-1))
      REAL*8 delt, const1
      REAL*8 mu00, Msolv
      REAL*8 cell_length, kv, ka, densityc, densitys
      REAL*8 r0, alpha, g, kg, b, kb
      REAL*8 moment0(NSTEP), moment1(NSTEP), moment2(NSTEP)
      REAL*8 moment3(NSTEP), moment4(NSTEP)
      REAL*8 time (NSTEP), concentration(NSTEP), seed_moment1(NSTEP)
      REAL*8 seed_moment2(NSTEP), seed_moment3(NSTEP)
      REAL*8 transmittance(NSTEP)
      REAL*8 temperature(NSTEP), relsatn(NSTEP)
      REAL*8 Temp, Csat, detFTVF
      REAL*8 derivtheta(NSTEP,NU,NTHETA)
      REAL*8 derivIg(NSTEP), derivconcg(NSTEP)
      REAL*8 derivIlnkg(NSTEP),derivconclnkg(NSTEP)
      REAL*8 derivIb(NSTEP), derivIlnkb(NSTEP)
      REAL*8 derivconcb(NSTEP), derivconclnkb(NSTEP)
      REAL*8 sigma_trans, sigma_conc
      REAL*8 AA(3,3), BB(3), gamma(3), lmin, lmax

*lsodes' parameters
      
      INTEGER itol, iopt, itask, istate, mf
      INTEGER lrw, liw, iwork(100)
      REAL*8  rtol, atol, rwork(1800)      



      EXTERNAL MOMENTS, MOMENTSJ, DIVPAG, DEVCSF, DLINRG

      COMMON /GROWTH_DATA/kg, g
      COMMON /BIRTH_DATA/kb, b 
      COMMON /EXP_DATA/r0, alpha, mu00
      COMMON Coeff

      DO 628 I = 1, Nvar
         Coeff(I)=X(I)
628   CONTINUE


*Simulation parameters
*     controller time step in minutes
      delt = 1.0D0
*     noise for transmittance measurement
      sigma_trans=0.009D0
*     noise for concentration measurement
      sigma_conc=0.0005D0

*Parameters for experimental set-up
*     cell length for spectrophotometer in millimeter
*     This was modified from that in Miller's thesis because
*     his value (2.0) did not agree with his simulation results
      cell_length=1.77D0
*     mass of solvent in grams, converted from 2000 gallons
      Msolv=7.57D6
*     volume shape factor (Appendix C in Miller)
      kv=1.0D0
*     area shape factor (Appendix C in Miller)
      ka=6.0D0
*     heat transfer coefficient multiplied by surface area
*     in calorie/minute/degree C 
*     density of solvent in g/cm^3 (solvent is water)
      densitys=1.0D0
*     density of crystal in g/cm^3 (Appendix C in Miller)
      densityc=2.11D0
*     seed size at nucleation
      r0=0.0D0
*     crystal density*volume shape factor,
*     in gram crystal/micron^3/particle
*     (alpha*L^3=mass of particle)
      alpha=kv*densityc*(1.0D-4)**3


*Initial conditions
*     initial concentration, g/g solvent
      concentration(1) = 0.493D0
*
*       The initial moments were computed using the following
*       population density function:
*
*       f_0(L)= aL^2 + b*L + c
*
*	The distribution function is assumed to be symmetrical
*	with the peak at L_bar.  The function is equal to zero 
*	at L=L_bar-w/2 and L=L_bar+w/2, where w is the width 
*	paramter. In the program THETA_(Ntemp1+1) is the total
*	seed, THETA(Ntemp1+2) is L_bar, and THETA_T(Ntemp1+3)
*	is the width.  Given the total mass, L_bar, and width w, 
*	the coefficient a, b, and c can be calcuted from the
*	following system of equlations.  Let
*	
*	lmin = L_bar - w/2
*	lmax = L_bar + w/2
*	mass = total seed mass 
*
*	Then
*
*	(lmax^6-lmin^6)/6 a + (lmax^5-lmin^5)/5 b + (lmax^4-lmin^4)/4 c =
*	mass/(mass_solvent*crystal_density)
*
*	lmax^2 a + lmax b + c = 0
*	lmin^2 a + lmax b + c = 0
*
*	Note: In the implementation, the first equation is scaled.
*

      lmin = x(Ntemp1+2)*1.0D3-
     &	     x(Ntemp1+3)*(1.0D-2)*x(Ntemp1+2)*1.0D3
      lmax = x(Ntemp1+2)*1.0D3+
     &       x(Ntemp1+3)*(1.0D-2)*x(Ntemp1+2)*1.0D3
*	print*,x(Ntemp1+1),x(Ntemp1+2),x(Ntemp1+3)
*	print*,lmin, lmax

      AA(1,1) = (lmax**6-lmin**6)/6.0D0/1.0D12 
      AA(1,2) = (lmax**5-lmin**5)/5.0D0/1.0D12
      AA(1,3) = (lmax**4-lmin**4)/4.0D0/1.0D12
      AA(2,1) = lmin**2
      AA(2,2) = lmin
      AA(2,3) = 1.0D0
      AA(3,1) = lmax**2
      AA(3,2) = lmax
      AA(3,3) = 1.0D0
      BB(1)=x(Ntemp1+1)*1.0D3/
     &		(Msolv*densityc)*(1.0D4)**3/1.0D12
      BB(2)=0.0D0
      BB(3)=0.0D0

*      print*,AA(1,1),AA(1,2),AA(1,3)
*      print*,AA(2,1),AA(2,2),AA(2,3)
*      print*,AA(3,1),AA(3,2),AA(3,3)

      CALL DLSARG (Norder, AA, LDA, BB, IPATH, gamma)
*      print*,gamma(1), gamma(2), gamma(3)



*     initial zeroth moment, number of particle/g solvent
      mu00 = gamma(1)*(lmax**3-lmin**3)/3.0D0+
     &	     gamma(2)*(lmax**2-lmin**2)/2.0D0+
     &	     gamma(3)*(lmax-lmin)
      moment0(1)= mu00
*     initial first moment, micron/g solvent
      moment1(1) = gamma(1)*(lmax**4-lmin**4)/4.0D0+
     &		   gamma(2)*(lmax**3-lmin**3)/3.0D0+
     &	           gamma(3)*(lmax**2-lmin**2)/2.0D0
      seed_moment1(1)=moment1(1)
*     initia1 second moment, micron^2/g solvent
      moment2(1) = gamma(1)*(lmax**5-lmin**5)/5.0D0+
     &		   gamma(2)*(lmax**4-lmin**4)/4.0D0+
     &	           gamma(3)*(lmax**3-lmin**3)/3.0D0
      seed_moment2(1)=moment2(1)
*     initial third moment, micron^3/g solvent
      moment3(1) = gamma(1)*(lmax**6-lmin**6)/6.0D0+
     &		   gamma(2)*(lmax**5-lmin**5)/5.0D0+
     &	           gamma(3)*(lmax**4-lmin**4)/4.0D0
      seed_moment3(1)=moment3(1)
*	print*,moment3(1)*Msolv*densityc*(1.0e-12)
*	print*,THETA_T(Ntemp1+1)
*     initial fourth moment, micron^4/g solvent
      moment4(1) = gamma(1)*(lmax**7-lmin**7)/7.0D0+
     &		   gamma(2)*(lmax**6-lmin**6)/6.0D0+
     &	           gamma(3)*(lmax**5-lmin**5)/5.0D0
*     initial relative supersaturation
      relsatn(1)=(concentration(1)-Csat(Temp(0.0D0)))/
     &     Csat(Temp(0.0D0))
*     initial transmittance measurement
      transmittance(1)=DEXP(-ka/2D0*cell_length/10.D0*moment2(1)*
     &     (densitys*(1.D-4)**2.D0))

*	print*,moment0(1)
*	print*,moment1(1)
*	print*,moment2(1)
*	print*,moment3(1)
*	print*,moment4(1)

*     initial d I/dg
      derivIg(1)=0.0D0
*     initial d conc/dg
      derivconcg(1) = 0.0D0
*     initial d I/ d lnkg
      derivIlnkg(1)=0.0D0
*     initial d conc/d lnkg
      derivconclnkg(1) = 0.0D0

*     initial d I/db
      derivIb(1)=0.0D0
*     initial d conc/db
      derivconcb(1) = 0.0D0
*     initial d I/ d lnkb
      derivIlnkb(1)=0.0D0
*     initial d conc/d lnkb
       derivconclnkb(1) = 0.0D0

*     Other initial conditions for derivatives wrt to theta
      DO 163 I = 1 , NU
         DO 164 J = 1 , NTHETA
            derivtheta(1,I,J)=0.0D0
164      CONTINUE
163   CONTINUE     



*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
      g=1.32D0
*     (mirons/minute)
      kg=DEXP(8.849D0)
*     (dimensionless)
      b=1.78D0
*     (number of particles/cm^3/minute) 
*     (the units have been corrected from that reported in
*     Table 3.1 in Miller)
      kb=DEXP(17.142D0)


*Simulation parameters
*********************************************************
      mf=222
      itask=1
      istate =1
      iopt=0
      lrw=3800
      liw=200
      rtol=1.0d-11
      atol=1.0d-10
      itol=1
***********************************************************

      time(1)=0.0D0
      T=0.0D0
      DO 1200 I=1,(NSTEP-1)
	 temperature(I)=Temp(T)
         Y(1)=moment0(I)
         Y(2)=moment1(I)
         Y(3)=moment2(I)
         Y(4)=moment3(I)
         Y(5)=moment4(I)
         Y(6)=concentration(I)
         Y(7)=seed_moment1(I)
         Y(8)=seed_moment2(I)
         Y(9)=seed_moment3(I)
         K = 10
         DO 263 M = 1 , NU
            DO 264 NNN = 1 , NTHETA
               Y(K)=derivtheta(I,M,NNN)
	       K = K + 1
264         CONTINUE
263      CONTINUE
*	Write(*,*)Coeff(1),Coeff(2)
        CALL lsodes ( MOMENTS,NEQ,Y,T,T+1,itol,rtol,
     &          atol,itask,istate,iopt, 
     &           rwork,lrw,iwork,liw,MOMENTSJ, mf )
         moment0(I+1)=Y(1)
         moment1(I+1)=Y(2)
         moment2(I+1)=Y(3)
         moment3(I+1)=Y(4)
         moment4(I+1)=Y(5)
         concentration(I+1)=Y(6)
         seed_moment1(I+1)=Y(7)
         seed_moment2(I+1)=Y(8)
         seed_moment3(I+1)=Y(9)
         relsatn(I+1)=(Y(6)-Csat(Temp(T)))/Csat(Temp(T))
         transmittance(I+1)=DEXP(-ka/2.0D0*cell_length/10.0D0*
     &        moment2(I+1)*(densitys*(1.0D-4)**2.D0))
         time(I+1)=T
         K = 10
         DO 363 M = 1 , NU
            DO 364 NNN = 1 , NTHETA
               derivtheta(I+1,M,NNN)=Y(K)
	       K = K + 1
364         CONTINUE
363      CONTINUE
         const1=ka/2.D0*cell_length/10.0D0*
     &			densitys*(1.0D-4)**(2.D0)