paramsun.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.
*
*      paramsun.f
*
*	This program reads the actual growth and nucleation 
*       kinetic parameters (g, kg, b, and kb), concentraion and 
*       transmittance data from a file, and estimates growth  
*       and nucleation kinetic parameters with the smallest
*       coefficient of variation
*      
*       The data files should be formatted in the following:
*
*                      actual_g         actual_Ln(kg)
*                      actual_b         actual_Ln(kb)
*                      concentration    transmittance
*
*       Ndata=Nset=the number of the data files
*
* 
*	The main program initializes the variables used by the
*	FFSQP subroutine.  Explanation of the variables is given 
*	in the FFSQP subroutine.  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, I
      PARAMETER(nparam = 4,nf=1, nineq=0, 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)
      REAL*8 bigbnd, eps, epseqn, udelta
      REAL*8 x(nparam), bl(nparam), bu(nparam)
      REAL*8 f(nf),g(nineq+neq+1),w(nwsize)
      INTEGER NSTEP,Nsets
      PARAMETER (NSTEP=161,Nsets=1)
      INTEGER NWRITE, flag

      REAL*8 ConcData(Nsets,NSTEP-1), TransData(Nsets,NSTEP-1)
*      REAL*8 actual_g, actual_ln_kg, actual_b, actual_ln_kb
      EXTERNAL FCN,cntr,grobfd,grcnfd

      COMMON /EXPERIMENAL_DATA/ConcData, TransData
      COMMON /DATA1/flag
      COMMON /DATA_OUT/NWRITE 

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

*  Initial guess:

       DATA x/	1.32D0, 8.84D0,
     &		1.78D0, 17.142D0/

* Lower bound:
  	bl(1)=1.0D0
        bl(2)=5.0D0
        bl(3)=1.0D0
        bl(4)=15.0D0
* Upper bound:
  	bu(1)=2.0D0
        bu(2)=10.0D0
        bu(3)=3.0D0
        bu(4)=20.0D0

	DO 121 I =1 , Nsets
      	   READ(55+I,*) actual_g, actual_ln_kg
           READ(55+I,*) actual_b, actual_ln_kb
* Read in data from file
          DO 122 J = 1,(NSTEP-1)
             READ(55+I,*) ConcData(I,J), TransData(I,J)
122       CONTINUE
          CLOSE(55+I)
121     CONTINUE




      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)


      STOP
      END

*****************************************************************

      SUBROUTINE cntr(Nvar,jjj,x,gj)
*
*	cntr is the subroutine for the constraints.  The
*	constraints are listed in the following order:   
*	Nonlinear inequality constraints, linear inequality 
*	constraints, nonlinear equality constraints, and linear 
*	equality constraints.
*
*     input:   Nvar - number of parameters
*	       jjj  - indicates the jjj_th constraint
*	       x    - Nvar-dimensional vector of parameters
*	       
*     output:  gj   - the jjj_th constraint
*

      INTEGER Nvar,jjj
      REAL*8 x(*),gj

      RETURN
      END


******************************************************************
*
*       SUBROUTINE FCN(Nvar,j,THETA,Phi)
*
*       Nvar = number of parameters to be determined
*       THETA = vector of length N, parameters to be determined
*       Phi = objective function to be optimized
*
*       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 purpose of this program is to serve as a benchmark
*       optimal control problem.  The objective of the optimal
*       control problem is to compute a temperature profile that
*       optimizes a crystal property of interest (see below).
*       The program reads the temperature profile from the
*       FORTRAN function Temp.  The output of the
*       program are plots of the moments, solute concentration,
*       relative supersaturation, temperature, and transmittance,
*       and values for 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.
*
*
*       Date:    February 12, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign
*
*
*

      SUBROUTINE FCN(Nvar,j,THETA,Phi)

      INTEGER NN, NEQ,j,Nvar,MXPARM,Ndata
      PARAMETER (NN=161, NEQ=9,MXPARM=50,Ndata=1)
      REAL*8 THETA(*), Phi
      INTEGER kfinal, I, Iteration
*      INTEGER INORM, METH, MITER
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
      REAL*8 T, Y(NEQ)
*      REAL*8 TOL, Hinit, Hmax, Hmin, MaxStep
      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 moment0(NN), moment1(NN), moment2(NN)
      REAL*8 moment3(NN), moment4(NN)
      REAL*8 time (NN), concentration(NN), seed_moment1(NN)
      REAL*8 seed_moment2(NN), seed_moment3(NN), transmittance(NN)
      REAL*8 temperature(NN), relsatn(NN)
      REAL*8 Temp, Csat
      REAL*8 ConcData(Ndata,NN-1), TransData(Ndata,NN-1)
      REAL*8 sigma_conc, sigma_trans
      REAL*8 weight_conc, weight_trans
      REAL*8 error_conc, error_trans
      REAL*8 DRNNOF
      REAL*8 AA(3,3), BB(3), gamma(3), lmin, lmax
      REAL*8 mass_seed, mean_seed, width
      
      
*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 /EXPERIMENAL_DATA/ConcData, TransData
      COMMON Iteration
 

      Phi=0.0D0
      error_conc=0.0D0
      error_trans=0.0D0
 
      DO 777 Iteration=1, Ndata
*Parameters for DIVPAG/DIVPRK
*     Tolerance
*      TOL = 1.0D0
*     Initial step size
*      Hinit = 1.0D-8
*     Maximum step size
*      Hmax = 1.50D-4
*     Minimum step size
*      Hmin = 0.0D0
*     Maximum number of steps	
*      MaxStep = 1.0D12
*     Error Norm
*      INORM = 1
*     Integration Method: 1 -Adams-Moulton; 2-Gear BDF
*      METH = 2 
*     Nonlinear Solveer: 0 - functional interation
*                        1 - chord method with user-defined
*                            Jacobian
*                        2 - chord method with divided -difference
*                            Jacobian
*       MITER = 2

*


*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

*     Total mass of seed crystals (grams)	
      mass_seed = 230.0D0
*     Mean size of the seed crystals (microns)
      mean_seed = 196.0D0
*     Percent width of the initial seed distribution
      width = 6.12D0



*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
      g=THETA(1)
*     (mirons/minute)
      kg=DEXP(THETA(2))
*     (dimensionless)
      b=THETA(3)
*     (number of particles/cm^3/minute) 
      kb=DEXP(THETA(4))


*Optimization parameters:
      sigma_trans=0.009D0
      sigma_conc=0.0005D0
      weight_trans=(1.0D0/sigma_trans)**2.0D0
      weight_conc=(1.0D0/sigma_conc)**2.0D0

*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 
*