📄 param2.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.
*
* param2.f
*
* This program reads measured concentration and moments data
* from files and estimates crystal growth and nucleation
* kinetic parameters using the weighted least-squares method.
* The data for this program are appropriate for the
* batch crystallization of crystals with one characteristic
* dimension.
*
* There are two input data files. The first file, named
* concen_data11, should contain the concentration data. In
* case you have second concentration data set, its name
* should be "concen_data12", and so on. The second file,
* named "mu_data11" should contain the measured moments data,
* in the format as follows:
*
* moment00, moment10, moment01, moment11, moment20
* moment01, moment21, moment12, moment22,
* moment30, moment31
*
* If you have a second set of moment data, then it should
* be named "mu_data12", and so on.
*
* If you are using this code for the ChE 391 course for a
* crystal with one or two characteristic dimensions, you
* should only change CSTEP, MSTEP, Nsets, Cinterval,
* Minterval, x, weights, and Temp:
*
* CSTEP = total points in concentration measurements
* MSTEP = total points in moments measurements
* Nsets = numbers of data sets
* Cinterval = time interval between two concentration
* measurements
* Minterval = time interval between two moment measurements
* x = initial guess of parameters
* weights = there are 6 weights to quantify the accuracy of
* the measurements. Each weight is set to the
* inverse of the measurement error variance.
* Normally it is very difficult to obtain an
* accurate estimate of the zeroth order moment
* and the fourth order moments, so the
* corresponding weights for these moments are
* normally set to zero.
* Temp = temperature profile (T(t))
*
* Parameter changes should be done both in the main
* program and in the subroutine FCN. Changes to Temp
* should be made in the subroutine Temp.
*
* If you are using this parameter estimation code for
* another process, then you would need to change the model
* equations as well as the above variables.
*
* Output of the program is the best-fit model parameters,
* which is the vector called "x" in the output screen. The
* other parameters on the output screen give information on
* the convergence of the parameter estimation algorithm,
* which can be ignored for most purposes. The code takes
* up to 10 minutes to run, depending on the speed of the
* computer and the quality of the initial parameter guesses.
*
* The main program initializes the variables used by the
* FFSQP subroutine, which minimizes the least-squared errors.
* 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.
* Our experience is that the program converges provided that
* you enter reasonable initial guesses.
*
* Date: February 12, 1998
* Authors: Serena H. Chung and Richard D. Braatz
* Department of Chemical Engineering
* University of Illinois at Urbana-Champaign
* Modified:March 4, 2000
* By: David L. Ma and Richard D. Braatz
*
use MSIMSL
INTEGER CSTEP, MSTEP, Nsets, Cinterval, Minterval
PARAMETER (CSTEP=160,Nsets=1, MSTEP=6)
PARAMETER (Cinterval=1,Minterval=30)
INTEGER nparam, npos1, npos2,I1
PARAMETER (nparam = 6, npos1=12, npos2=8)
CHARACTER*60 file_concen, file_mu, bob
INTEGER nf, nineqn, nineq, neqn, neq, iwsize, nwsize
INTEGER mode, iprint,miter
INTEGER inform, I
PARAMETER(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 NWRITE, flag
REAL*8 ConcData(Nsets,CSTEP), Mu00Data(Nsets,MSTEP)
REAL*8 Mu10Data(Nsets,MSTEP), Mu01Data(Nsets,MSTEP)
REAL*8 Mu11Data(Nsets,MSTEP), Mu20Data(Nsets,MSTEP)
REAL*8 Mu02Data(Nsets,MSTEP), Mu21Data(Nsets,MSTEP)
REAL*8 Mu12Data(Nsets,MSTEP), Mu22Data(Nsets,MSTEP)
REAL*8 Mu30Data(Nsets,MSTEP), Mu31Data(Nsets,MSTEP)
REAL*8 weight_conc, weight_mu00, weight_mu10, weight_mu01,
& weight_mu11, weight_mu20, weight_mu02,
& weight_mu21, weight_mu12, weight_mu22,
& weight_mu30, weight_mu31
EXTERNAL FCN,cntr,grobfd,grcnfd
COMMON /CONCEN/ConcData
COMMON /MU00/Mu00Data
COMMON /MU10/Mu10Data
COMMON /MU01/Mu01Data
COMMON /MU11/Mu11Data
COMMON /MU20/Mu20Data
COMMON /MU02/Mu02Data
COMMON /MU21/Mu21Data
COMMON /MU12/Mu12Data
COMMON /MU22/Mu22Data
COMMON /MU30/Mu30Data
COMMON /MU31/Mu31Data
COMMON /DATA1/flag
COMMON /DATA_OUT/NWRITE
COMMON /WEIGHT/weight_conc,
& weight_mu00, weight_mu10, weight_mu01,
& weight_mu11, weight_mu20, weight_mu02,
& weight_mu21, weight_mu12, weight_mu22,
& weight_mu30, weight_mu31
* Initial guess:
DATA x/ 1.15D0, 7.5D0, 1.65D0, 8.8D0,
& 1.8D0, 18.174D0/
print*, x(1),x(2),x(3),x(4),x(5),x(6)
* file name
file_concen='concen_data'
file_mu='mu_data'
*Input the weights for least-squares parameter estimation
weight_conc=1.0D0/0.49D0/0.49D0
weight_mu00=0.0D0
weight_mu10=1.0D0/59008.0D0/59008.0D0
weight_mu01=1.0D0/199721.0D0/199721.0D0
weight_mu11=1.0D0/44892500.0D0/44892500.0D0
weight_mu20=1.0D0/15474300.0D0/15474300.0D0
weight_mu02=1.0D0/137020000.0D0/137020000.0D0
weight_mu21=0.0D0
weight_mu12=0.0D0
weight_mu22=0.0D0
weight_mu30=0.0D0
weight_mu31=0.0D0
*FFSQP parameters
bigbnd=1.0D12
eps=1.0D-8
epseqn=0.0D0
udelta=0.0D0
iprint=2
nineqn=0
neqn=0
* Lower bounds:
bl(1)=1.0D0
bl(2)=5.0D0
bl(3)=1.0D0
bl(4)=5.0D0
bl(5)=1.0D0
bl(6)=3.0D0
* Upper bounds:
bu(1)=2.0D0
bu(2)=15.0D0
bu(3)=2.0D0
bu(4)=15.0D0
bu(5)=3.0D0
bu(6)=20.0D0
DO I =1+10 , Nsets+10
I1=I-10
* Read in data from file
DO J = 1,CSTEP
WRITE(file_concen(npos1:npos1+1), '(I2)') I
file_concen = file_concen(1:npos1+1) // '.txt'
OPEN(UNIT=20, FILE=file_concen, FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='OLD')
READ(20,*) ConcData(I1,J)
ENDDO
CLOSE(UNIT=20)
DO J=1,MSTEP
WRITE(file_mu(npos2:npos2+1), '(I2)') I
file_mu = file_mu(1:npos2+1) // '.txt'
OPEN(UNIT=20, FILE=file_mu, FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='OLD')
READ(20,*) Mu00Data(I1,J), Mu10Data(I1,J), Mu01Data(I1,J),
& Mu11Data(I1,J), Mu20Data(I1,J), Mu02Data(I1,J),
& Mu21Data(I1,J), Mu12Data(I1,J), Mu22Data(I1,J),
& Mu30Data(I1,J), Mu31Data(I1,J)
ENDDO
CLOSE(UNIT=20)
ENDDO
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='param.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
WRITE(20,10)x(1)
WRITE(20,10)x(2)
WRITE(20,10)x(3)
WRITE(20,10)x(4)
WRITE(20,10)x(5)
WRITE(20,10)x(6)
10 FORMAT(E24.16)
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.
*
*
* Date: February 12, 1998
* Authors: Serena H. Chung and Richard D. Braatz
* Department of Chemical Engineering
* University of Illinois at Urbana-Champaign
* Modified:March 4, 2000
* By: David L. Ma and Richard D. Braatz
*
SUBROUTINE FCN(Nvar,j,THETA,Phi)
INTEGER CSTEP, MSTEP, Nsets, Cinterval, Minterval
PARAMETER (CSTEP=160,Nsets=1,MSTEP=6)
PARAMETER (Cinterval=1,Minterval=30)
INTEGER NN, NEQ,j,Nvar,MXPARM,Ndata
PARAMETER (NN=CSTEP+1,NEQ=19,MXPARM=50,Ndata=Nsets)
REAL*8 THETA(*), Phi
INTEGER kfinal, I, Iteration
INTEGER Norder, LDA, LDB, IPATH
PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
REAL*8 T, Y(NEQ)
REAL*8 delt, tfinal
REAL*8 mu00
REAL*8 cell_length, Msolv, kv, ka, UA, densityc, densitys
REAL*8 r0, alpha, g1, g2, kg1, kg2, b, kb
REAL*8 moment00(NN), moment10(NN), moment01(NN)
REAL*8 moment11(NN), moment20(NN), moment02(NN)
REAL*8 moment21(NN), moment12(NN), moment22(NN)
REAL*8 moment30(NN), moment31(NN)
REAL*8 time (NN), concentration(NN)
REAL*8 concentration_measured(NN)
REAL*8 seed_moment10(NN), seed_moment01(NN)
REAL*8 seed_moment11(NN), seed_moment20(NN)
REAL*8 seed_moment02(NN), seed_moment21(NN)
REAL*8 seed_moment12(NN)
REAL*8 temperature(NN), relsatn(NN)
REAL*8 Temp, Csat
REAL*8 ConcData(Ndata,CSTEP), Mu00Data(Nsets,MSTEP)
REAL*8 Mu10Data(Nsets,MSTEP), Mu01Data(Nsets,MSTEP)
REAL*8 Mu11Data(Nsets,MSTEP), Mu20Data(Nsets,MSTEP)
REAL*8 Mu02Data(Nsets,MSTEP), Mu12Data(Nsets,MSTEP)
REAL*8 Mu21Data(Nsets,MSTEP), Mu22Data(Nsets,MSTEP)
REAL*8 Mu30Data(Nsets,MSTEP), Mu31Data(Nsets,MSTEP)
REAL*8 weight_conc, weight_mu00, weight_mu10
REAL*8 weight_mu01, weight_mu11, weight_mu20
REAL*8 weight_mu02, weight_mu21, weight_mu12
REAL*8 weight_mu22, weight_mu30, weight_mu31
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, DLSARG
COMMON /GROWTH_DATA1/kg1, g1
COMMON /GROWTH_DATA2/kg2, g2
COMMON /BIRTH_DATA/kb, b
COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv
COMMON Iteration
COMMON /CONCEN/ConcData
COMMON /MU00/Mu00Data
COMMON /MU10/Mu10Data
COMMON /MU01/Mu01Data
COMMON /MU11/Mu11Data
COMMON /MU20/Mu20Data
COMMON /MU02/Mu02Data
COMMON /MU21/Mu21Data
COMMON /MU12/Mu12Data
COMMON /MU22/Mu22Data
COMMON /MU30/Mu30Data
COMMON /MU31/Mu31Data
COMMON /WEIGHT/weight_conc,
& weight_mu00, weight_mu10, weight_mu01,
& weight_mu11, weight_mu20, weight_mu02,
& weight_mu21, weight_mu12, weight_mu22,
& weight_mu30, weight_mu31
Phi=0.0D0
error_conc=0.0D0
error_trans=0.0D0
DO 777 Iteration=1, Ndata
*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)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -