# [depends] flow.dat lc.dat# This py-file loads data file 'flow.dat'.# This py-file loads data file 'lc.dat'.fromparestoimportParestofrommiscimportgnuplotsaveimportnumpyasnpfromscipy.interpolateimportinterp1dimportmatplotlib.pyplotaspltimportcasadiimportosplt.ion()# Initialize modelpe=Paresto()pe.print_level=1pe.transcription='shooting'pe.x=['VR','nA','nB','nC','nD']pe.p=['k1','k2','cBf']pe.nlp_solver_options={'ipopt':{'mumps_scaling':0}}# Dependent variables with definitionspe.y=['lc']pe.h=lambdat,v,p:[1/(1+2*v.nD/casadi.fmax(v.nC,1e-6))]# avoid divide-by-zero# Data and measurementspe.d=['Qf','lc_m']# ODE right-hand-sidepe.ode=lambdat,v,p:[v.Qf,-p.k1*v.nA*v.nB/v.VR,v.Qf*p.cBf-v.nB*(p.k1*v.nA+p.k2*v.nC)/v.VR,v.nB*(p.k1*v.nA-p.k2*v.nC)/v.VR,p.k2*v.nC*v.nB/v.VR]# Relative least squares objective functionpe.lsq=lambdat,y,p:[y.lc_m/y.lc-1]# Load datateaf=0.00721teaden=0.728flow=np.loadtxt('flow.dat')lc=np.loadtxt('lc.dat')tQf=np.concatenate(([0.],flow[:,0]))Qf=np.concatenate(([0.],flow[:,1]/teaden))tlc=lc[:,0]lc=lc[:,1]# Get all time points occurring in either tlc or tQftout,ic=np.unique(np.concatenate((tQf,tlc)),return_inverse=True)Qf_ind=ic[:len(tQf)]lc_ind=ic[len(tQf):]# Interpolate lcmeas and Qf to this new gridQf=interp1d(tQf,Qf,kind='previous',fill_value="extrapolate")(tout)lc_m=interp1d(tlc,lc,kind='previous',fill_value="extrapolate")(tout)N=Qf.shape[0]# Replace NaNs with zerosQf[np.isnan(Qf)]=0.lc_m[np.isnan(lc_m)]=0.# Initial volumeVR0=2370# Optionspe.tout=toutpe.lsq_ind=lc_ind# only include lc_ind in objective# must tell paresto number of meas points since pe.d# has 2 (flow, meas) instead of 1 (meas) elementspe.ndata=len(pe.lsq_ind)# Create a paresto instancepe.initialize()# Solution in the booknA0=2.35k1=2500k2=1250# Initial guess for parameterspe.ic.k1=k1pe.ic.k2=0.9*k2pe.ic.cBf=teafpe.ic.VR=VR0pe.ic.nA=1.1*nA0pe.ic.nB=0pe.ic.nC=0pe.ic.nD=0pe.lb_ub()pe.lb.k1=0.5*pe.ic.k1pe.lb.k2=0.5*pe.ic.k2pe.lb.nA=0.5*nA0pe.ub.k1=1.5*pe.ic.k1pe.ub.k2=1.5*pe.ic.k2pe.ub.nA=1.5*nA0# Estimate parametersest,v,p=pe.optimize(np.vstack([Qf,lc_m]).reshape(2,N,1))# Also calculate confidence intervals with 95 % confidenceconf=pe.confidence(est,0.95)print('Estimated parameters')print('\n'.join(f'{k}: {float(v)}'fork,vinest.theta.items()))print('Bounding box intervals')print('\n'.join(f'{k}: {float(v)}'fork,vinconf.bbox.items()))# Save data (matches Octave: data.model = [tout, states..., lc], data.measurement = [tlc, lc])model_cols=[pe.tout]+[v[s].flatten()forsinpe.x]+[v['lc'].flatten()]data={'model':np.column_stack(model_cols),'measurement':np.column_stack((tlc,lc)),}gnuplotsave('bvsm.dat',data)# Write code for saving dataifnotos.getenv('OMIT_PLOTS')=='true':plt.figure(figsize=(10,8))plt.subplot(2,2,1)plt.plot(pe.tout,v['nA'].flatten(),label='n_A')plt.plot(pe.tout,v['nC'].flatten(),label='n_C')plt.plot(pe.tout,v['nD'].flatten(),label='n_D')plt.legend()plt.xlabel('time (min)')plt.ylabel('Amount of substance (kmol)')plt.title('Amount of substance of species A, C, and D versus time')plt.grid(True)plt.subplot(2,2,2)plt.plot(pe.tout,v['nB'].flatten(),label='n_B')plt.legend()plt.xlabel('time (min)')plt.ylabel('Amount of substance (kmol)')plt.title('Amount of substance of species B versus time')plt.grid(True)plt.subplot(2,2,3)plt.step(pe.tout,v['Qf'].flatten(),where='mid')plt.xlabel('time (min)')plt.ylabel('flowrate (kg/min)')plt.title('Base addition rate')plt.grid(True)plt.subplot(2,2,4)plt.plot(pe.tout,v['lc'].flatten(),label='model')plt.plot(tlc,lc,'o',label='measurement')plt.ylim(0,2*max(lc))plt.legend()plt.xlabel('time (min)')plt.title('LC measurement')plt.grid(True)plt.tight_layout()plt.show()
