📄 optexp.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.
*
* optexp.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 (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 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+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),g1(nineq+neq),w(nwsize)
REAL*8 g,kg,b,kb
REAL*8 scale(nparam)
EXTERNAL FCN,cntr,grobfd,grcnfd
COMMON /GROWTH_DATA/kg, g
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 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
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) = 95.0D0
* Read g, kg, b, kb
OPEN(UNIT=20, FILE='param.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='OLD')
READ(20,*)g
READ(20,*)kg
READ(20,*)b
READ(20,*)kb
CLOSE(UNIT=20)
kg=DEXP(kg)
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.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)
* 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)=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,g1,
* 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+3)
INTEGER NSTEP, NEQ, MXPARM, NTHETA, NU
PARAMETER (NSTEP=161,NEQ=33,MXPARM=50,NTHETA=4, NU=6)
INTEGER I, J, K, M, N
INTEGER Norder, LDA, LDB, IPATH
PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
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 delt, tfinal, Msolv
REAL*8 mu00
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 temperature(NSTEP), relsatn(NSTEP)
REAL*8 Temp, Csat, detFTVF
REAL*8 derivtheta(NSTEP,NU,NTHETA)
REAL*8 conc_variance, mu0_variance, mu1_variance
REAL*8 mu3_variance, mu4_variance, mu2_variance
REAL*8 AA(3,3), BB(3), gamma1(3), lmin, lmax
*lsodes' parameters
INTEGER itol, iopt, itask, istate, mf
INTEGER lrw, liw, iwork(900)
REAL*8 rtol, atol, rwork(5000)
EXTERNAL MOMENTS, MOMENTSJ
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
* print*, "in Fcn"
*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)
* 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)
*Noise estimates
mu0_variance=0.1456D0
mu1_variance=9.3523D0
mu2_variance=3484.3446D0
mu3_variance=1.6221D6
mu4_variance=7.043D8
conc_variance=3.53D-4
*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
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
CALL DLSARG (Norder, AA, LDA, BB, IPATH, gamma1)
* initial zeroth moment, number of particle/g solvent
mu00 = gamma1(1)*(lmax**3-lmin**3)/3.0D0+
& gamma1(2)*(lmax**2-lmin**2)/2.0D0+
& gamma1(3)*(lmax-lmin)
moment0(1)= mu00
* initial first moment, micron/g solvent
moment1(1) = gamma1(1)*(lmax**4-lmin**4)/4.0D0+
& gamma1(2)*(lmax**3-lmin**3)/3.0D0+
& gamma1(3)*(lmax**2-lmin**2)/2.0D0
seed_moment1(1)=moment1(1)
* initia1 second moment, micron^2/g solvent
moment2(1) = gamma1(1)*(lmax**5-lmin**5)/5.0D0+
& gamma1(2)*(lmax**4-lmin**4)/4.0D0+
& gamma1(3)*(lmax**3-lmin**3)/3.0D0
seed_moment2(1)=moment2(1)
* initial third moment, micron^3/g solvent
moment3(1) = gamma1(1)*(lmax**6-lmin**6)/6.0D0+
& gamma1(2)*(lmax**5-lmin**5)/5.0D0+
& gamma1(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) = gamma1(1)*(lmax**7-lmin**7)/7.0D0+
& gamma1(2)*(lmax**6-lmin**6)/6.0D0+
& gamma1(3)*(lmax**5-lmin**5)/5.0D0
* initial relative supersaturation
relsatn(1)=(concentration(1)-Csat(Temp(0.0D0)))/
& Csat(Temp(0.0D0))
print*, moment0(1)
* 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
*Simulation parameters
*********************************************************
*
mf=222
itask=1
istate =1
iopt=0
lrw=5000
liw=900
rtol=1.0d-11
atol=1.0d-10
itol=1
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -