optconsun.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.
*
*	optcon.f
*	
*
*	This program calculates the optimum temperature profile
*	and initial seed distribution for optimizing a final time 
*	crystal	property of interest.  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:    June 19, 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=2, 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-6
      epseqn=0.0D0
      udelta=0.0D0
      iprint=2
      nineqn=1
      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.5166629898298D-02
      x(2) = -2.0218060663727D-02
      x(3) = -2.9765463427082D-02
      x(4) = -4.5824907459300D-02
      x(5) = -8.7417830856041D-02
      x(6) = -1.0000000000000D-01
      x(7) = -1.0000000000000D-01
      x(8) = -1.0000000000000D-01
      x(9) =  20.0D0
      x(10) = 5.0D0
      x(11) = 5.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)


*      DATA x/	-0.24679704220759D-01, -0.24679704220759D-01,
*     &		-0.17028909974573D-33, -0.17028909974573D-33, 
*     &		-0.26141823689473D-01, -0.26141823689473D-01,
*     &		-0.48379245469584D-01, -0.48379245469584D-01,
*     &		-0.10000000000000D+00, -0.10000000000000D+00,
*     &		-0.10000000000000D+00, -0.10000000000000D+00,
*     &		-0.10000000000000D+00, -0.10000000000000D+00,
*     &		-0.10000000000000D+00, -0.10000000000000D+00/

* 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(Nvar,j,THETA_T,function)
*
*	Nvar		number of parameters
*	j		indicates the jth function (not 
*			relevent to the current optimization
*			probem)
*	THETA_T		Nvar-dimensional vector of parameters
*       function	objective function value
*	
*
*       This subroutine is a version of the program open.f. 
*	The 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 objective function can be any of the three common final
*	crystal properties of interest: (i) weight mean size, 
*	(ii) ratio of nucleated crystal mass to seed crystal mass, 
*	and (iii) coefficient of variance.
*
*	The following subroutines are required for FCN:
*	
*		Temp
*		MOMENTS
*		MOMENTSJ
*		Csat
*		growth
*		birth
*

      SUBROUTINE FCN(Nvar,j,THETA_T,function)

      INTEGER NN, NEQ, MXPARM,Nvar,j, Ntemp1
      PARAMETER (NN=161, NEQ=9, MXPARM=50,Ntemp1 = 8)
      INTEGER  kfinal, I
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
      REAL*8 T, Y(NEQ)
      REAL*8 THETA_T(*), function, Coeff(Ntemp1+3)
      REAL*8 delt, tfinal
      REAL*8 mu00
      REAL*8 cell_length, Msolv, kv, ka, UA, densityc, densitys
      REAL*8 r0, alpha, g, kg, b, kb
      REAL*8 weight_mean_size, cov, mass_ratio
      REAL*8 moment0(NN), moment1(NN), moment2(NN)
      REAL*8 moment3(NN), moment4(NN)
      REAL*8 time (NN), concentration(NN), seed_moment1(NN)
      REAL*8 concentration_measured(NN)
      REAL*8 seed_moment2(NN), seed_moment3(NN), transmittance(NN)
      REAL*8 temperature(NN), relsatn(NN)
      REAL*8 Temp, Csat
      REAL*8 AA(3,3), BB(3), gamma(3), lmin, lmax
      
      
*lsodes' parameters
      
      INTEGER itol, iopt, itask, istate, mf
      INTEGER lrw, liw, iwork(200)
      REAL*8  rtol, atol, rwork(3800)      
      
      
      
	
      EXTERNAL MOMENTS, MOMENTSJ, DIVPAG,DRNNOF, DLSARG
      COMMON /GROWTH_DATA/kg, g
      COMMON /BIRTH_DATA/kb, b 
      COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv
      COMMON Coeff
      COMMON concentration

      DO 1111 I = 1, Nvar
      		Coeff(I)=THETA_T(I)
1111  CONTINUE

*Simulation parameters
*     controller time step in minutes
      delt = 1.0D0
*     total time step
      kfinal=NN
*     final time in minutes
      tfinal = DFLOAT(NN)*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



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

*
*       The initial moments were computed using the following
*       population density function:
*
*       f0(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.
*	Note: L_bar = (first moment)/(zeroth moment)
*

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

      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)=THETA_T(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)