optexp2.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.
*
*	optexp2.f
*	
*
*	This program calculates the temperature profile and
*	initial seed distribution for optimizing the 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:
*
*		Subroutine/Program	Parameter
*		------------------	---------
*		Main			ntemp
*		FCN			Ntemp1
*		cntr			Ntemp2
*		Temp			Ntemp3 
*
*	The initial seed distribution is characterized by the
*	total mass, mean size, and width of the distribution.
*
*	Parameter inputs are "best guesses" for the growth and 
*       nucleation kinetic parameters (g1, kg1, g2, kg2, 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 included.
*	
*
*       Date:    July 6, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*	Modified:March 20, 2000
*	By:      David L. Ma and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign
*
************************************************************************
*      PROGRAM MAIN
	 use MSIMSL

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

      INTEGER nparam, nf, nineqn, nineq, neqn, neq, iwsize, nwsize
      INTEGER mode, iprint,miter
      INTEGER inform, ntemp
      PARAMETER(ntemp = 8, nparam=ntemp+4,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 g1,kg1,g2,kg2,b,kb
      REAL*8 scale(nparam)
      EXTERNAL FCN,cntr,grobfd,grcnfd
	COMMON /GROWTH_DATA1/kg1, g1
	COMMON /GROWTH_DATA2/kg2, g2
      COMMON /BIRTH_DATA/kb, b

      bigbnd=1.0D12
      eps=1.0D-4
      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 moment10 of the seed distribution 
* in micron/g solvent.  x(ntemp+3) is moment01.  
* x(ntemp+4) 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) = 600.0D0
*	 x(12) = 95.0D0

      x(1) =-0.1D0
      x(2) =-0.1D0
      x(3) =-0.1D0
      x(4) =-0.1D0
      x(5) =-0.1D-8
      x(6) =-0.1D-8
      x(7) =-0.1D0
      x(8) =-0.1D-8
      x(9) = 5.0D0
      x(10) = 600.0D0
      x(11) = 600.0D0
	x(12) = 95.0D0

*	Read g1, kg1, g2, kg2, b, kb  
	OPEN(UNIT=20, FILE='param.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='OLD')
	READ(20,*)g1
	READ(20,*)kg1
	READ(20,*)g2
	READ(20,*)kg2 
      READ(20,*)b
	READ(20,*)kb
	CLOSE(UNIT=20)
	kg1=DEXP(kg1)
	kg2=DEXP(kg2)
      kb=DEXP(kb)
	

	
*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.0D-3
	  scale(ntemp+4)=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.5D0
*	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)
	bl(ntemp+4)=5.0D0*scale(ntemp+4)

* Upper bound for the parameters:
      DO 667 I = 1, ntemp
667	bu(I)=-1.0D-9
      bu(ntemp+1)=110000.0D0*scale(ntemp+1)
      bu(ntemp+2)=600.0D0*scale(ntemp+2)
      bu(ntemp+3)=600.0D0*scale(ntemp+3)
	bu(ntemp+4)=95.0D0*scale(ntemp+4)
      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)

	OPEN(UNIT=20, FILE='slope.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      PRINT*,'Final Solution'
      DO I = 1, nparam 
      	PRINT*,x(I)/scale(I)
	    WRITE(20,700)x(I)/scale(I)
	ENDDO
700	FORMAT(E27.16)
      PRINT*,'objective = ', f


	CLOSE(UNIT=20)
	

      STOP
      END

******************************************************************
      SUBROUTINE FCN(Nvar,jjj,x,fj)
      INTEGER Nvar,jjj,Ntemp1	
      PARAMETER (Ntemp1=8)
      REAL*8 x(*), fj
      REAL*8 Coeff(Ntemp1+4)
      INTEGER NSTEP, NEQ, MXPARM, NTHETA, NU
      PARAMETER (NSTEP=161,NEQ=91,MXPARM=50,NTHETA=6, NU=12)
      INTEGER I, J, K, M, N, P
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
      PARAMETER(LDFAC=6, P=6)
	INTEGER IPVT(P)
      REAL*8 T, Y(NEQ)
      REAL*8 F(NU*(NSTEP-1), NTHETA)
      REAL*8 FTVF(NTHETA, NTHETA)
      REAL*8 FTV(NTHETA,NU*(NSTEP-1))
      REAL*8 DET1, DET2, FAC(LDFAC,LDFAC)
	
      REAL*8 delt, tfinal, Msolv
      REAL*8 mu00
      REAL*8 cell_length, kv, ka, densityc, densitys
      REAL*8 r0, alpha, g1, kg1, g2, kg2, b, kb
      
	REAL*8 moment00(NSTEP), moment10(NSTEP), moment01(NSTEP)
      REAL*8 moment11(NSTEP), moment20(NSTEP), moment02(NSTEP)
      REAL*8 moment21(NSTEP), moment12(NSTEP), moment22(NSTEP)
      REAL*8 moment30(NSTEP), moment31(NSTEP)

      REAL*8 time (NSTEP), concentration(NSTEP) 
      REAL*8 seed_moment10(NSTEP), seed_moment01(NSTEP)
      REAL*8 seed_moment11(NSTEP), seed_moment20(NSTEP)
      REAL*8 seed_moment02(NSTEP), seed_moment21(NSTEP) 
      REAL*8 seed_moment12(NSTEP)

      REAL*8 temperature(NSTEP), relsatn(NSTEP)
      REAL*8 Temp, Csat, detFTVF
      REAL*8 derivtheta(NSTEP,NU,NTHETA)
      
	REAL*8 conc_variance, mu00_variance, mu10_variance
      REAL*8 mu01_variance, mu11_variance, mu20_variance
      REAL*8 mu02_variance, mu21_variance, mu12_variance
	REAL*8 mu22_variance, mu30_variance, mu31_variance
	

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

      EXTERNAL MOMENTS, MOMENTSJ
      COMMON /GROWTH_DATA1/kg1, g1
	COMMON /GROWTH_DATA2/kg2, g2
      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
*     final time in minutes
      tfinal = DFLOAT(NSTEP)*delt

*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

*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
*      g1=1.32D0
*     (mirons/minute)
*      kg1=DEXP(8.849D0)
*     (dimensionaless)
*      g2=1.32D0
*     (mirons/minute)
*      kg2=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)



*Initial conditions
*     initial concentration, g/g solvent
      concentration(1) = 0.493D0


*       The initial moments were computed using the following
*       population density function:
*
*       f0= -0.00034786x^2 + 0.1363609*x - 13.2743
*	    -0.00034786y^2 + 0.1363609*y - 13.2743
*
*       which is in units of number of crystals/g solvent/micron.
*       This function was based on assuming 0.05 g as the mass in
*       1650 grams of solvent and 180-212 microns as the size range
*       for the seed crystals, and assuming a quadratic distibution
*       function.  The values are valid if mass of seed crystals
*       is scaled up proportionally to the mass of solvent.
*
*     initial zeroth moment, number of particle/g solvent
      mu00=121.575D0
      moment00(1)=mu00
           
*     initial moment10, moment01 micron/g solvent
      moment10(1)=2.383D4
      moment01(1)=moment10(1)
      seed_moment10(1)=moment10(1)
      seed_moment01(1)=moment10(1)
      
*     initia1 moment11, micron^2/g solvent
      moment11(1)=4.67D6
      seed_moment11(1)=moment11(1)
      
*     initial moment20, moment02 micron^2/g solvent
      moment20(1)=4.679D6
      moment02(1)=moment20(1)
      seed_moment20(1)=moment20(1)
      seed_moment02(1)=moment20(1)
      
*     initial moment12, moment21 micron^3/g solvent
      moment12(1)=9.17D8
      moment21(1)=moment12(1)
      seed_moment21(1)=moment12(1)
      seed_moment12(1)=moment12(1)