From a7d9785040b127e4f683077344f70a32aee79bc5 Mon Sep 17 00:00:00 2001 From: schenkch Date: Mon, 3 Jun 2019 17:30:05 -0400 Subject: [PATCH] Updated example 3 --- .../paper_example3_generatesimabsdata.py | 62 +++++++++---------- .../paper_example3_inputsandnonabs.py | 57 +++++++++-------- 2 files changed, 59 insertions(+), 60 deletions(-) diff --git a/kipet/examples/PaperExamples/paper_example3_generatesimabsdata.py b/kipet/examples/PaperExamples/paper_example3_generatesimabsdata.py index 2680adbb..aeeac72d 100644 --- a/kipet/examples/PaperExamples/paper_example3_generatesimabsdata.py +++ b/kipet/examples/PaperExamples/paper_example3_generatesimabsdata.py @@ -83,7 +83,7 @@ def Lorentzian_parameters(): params_g['gammas'] = [100.0, 30000.0, 100.0, 100.0] - return {'AH': params_d,'B':params_a,'C': params_c, 'BH+': params_b, 'A-': params_f, 'P':params_e, 'AC-': params_g} + return {'A': params_d,'B':params_a,'C': params_c, 'D': params_b, 'E': params_f, 'G':params_e, 'F': params_g} def frange(start, stop, step): return takewhile(lambda x: x< stop, count(start, step)) @@ -108,13 +108,13 @@ def frange(start, stop, step): # components components = dict() - components['AH'] = 0.395555 + components['A'] = 0.395555 components['B'] = 0.0351202 components['C'] = 0.0 - components['BH+'] = 0.0 - components['A-'] = 0.0 - components['AC-'] = 0.0 - components['P'] = 0.0 + components['D'] = 0.0 + components['E'] = 0.0 + components['F'] = 0.0 + components['G'] = 0.0 builder.add_mixture_component(components) @@ -126,7 +126,7 @@ def frange(start, stop, step): params = dict() - params['k0'] = 49.7796 + params['k0'] = 0.2545 params['k1'] = 8.93156 params['k2'] = 1.31765 params['k3'] = 0.310870 @@ -142,21 +142,21 @@ def frange(start, stop, step): # stoichiometric coefficients gammas = dict() - gammas['AH'] = [-1, 0, 0, -1, 0] + gammas['A'] = [-1, 0, 0, -1, 0] gammas['B'] = [-1, 0, 0, 0, 1] gammas['C'] = [0, -1, 1, 0, 0] - gammas['BH+'] = [1, 0, 0, 0, -1] - gammas['A-'] = [1, -1, 1, 1, 0] - gammas['AC-'] = [0, 1, -1, -1, -1] - gammas['P'] = [0, 0, 0, 1, 1] + gammas['D'] = [1, 0, 0, 0, -1] + gammas['E'] = [1, -1, 1, 1, 0] + gammas['F'] = [0, 1, -1, -1, -1] + gammas['G'] = [0, 0, 0, 1, 1] def rule_algebraics(m, t): r = list() - r.append(m.Y[t, '0'] - m.P['k0'] * m.Z[t, 'AH'] * m.Z[t, 'B']) - r.append(m.Y[t, '1'] - m.P['k1'] * m.Z[t, 'A-'] * m.Z[t, 'C']) - r.append(m.Y[t, '2'] - m.P['k2'] * m.Z[t, 'AC-']) - r.append(m.Y[t, '3'] - m.P['k3'] * m.Z[t, 'AC-'] * m.Z[t, 'AH']) - r.append(m.Y[t, '4'] - m.P['k4'] * m.Z[t, 'AC-'] * m.Z[t, 'BH+']) + r.append(m.Y[t, '0'] - m.P['k0'] * m.Z[t, 'A'] * m.Z[t, 'B']) + r.append(m.Y[t, '1'] - m.P['k1'] * m.Z[t, 'E'] * m.Z[t, 'C']) + r.append(m.Y[t, '2'] - m.P['k2'] * m.Z[t, 'F']) + r.append(m.Y[t, '3'] - m.P['k3'] * m.Z[t, 'F'] * m.Z[t, 'A']) + r.append(m.Y[t, '4'] - m.P['k4'] * m.Z[t, 'F'] * m.Z[t, 'D']) return r #: there is no ae for Y[t,5] because step equn under rule_odes functions as the switch for the "C" equation @@ -185,10 +185,10 @@ def rule_odes(m, t): feed_times=[101.035, 303.126]#, 400. builder.add_feed_times(feed_times) - model = builder.create_pyomo_model(0, 600) + model = builder.create_pyomo_model(0, 600.) builder.add_absorption_data(S_frame) - write_absorption_data_to_txt('Sij_FEcaseexample5.txt', S_frame) + write_absorption_data_to_txt('Sij_FEcaseexample5-2.txt', S_frame) model = builder.create_pyomo_model(0., 600.) @@ -223,10 +223,10 @@ def rule_odes(m, t): #to this function as an argument dictionary #New Inputs for discrete feeds - Z_step = {'AH': .03} #Which component and which amount is added + Z_step = {'A': .03} #Which component and which amount is added X_step = {'V': .01} jump_states = {'Z': Z_step, 'X': X_step} - jump_points1 = {'AH': 101.035}#Which component is added at which point in time + jump_points1 = {'A': 101.035}#Which component is added at which point in time jump_points2 = {'V': 303.126} jump_times = {'Z': jump_points1, 'X': jump_points2} @@ -242,13 +242,13 @@ def rule_odes(m, t): #results_sim = sim.run_sim('ipopt', #tee=True, #solver_opts=options) - sigmas={'AH': 1e-10, + sigmas={'A': 1e-10, 'B':1e-10, 'C': 1e-10, - 'BH+': 1e-10, - 'A-':1e-10, - 'AC-':1e-10, - 'P':1e-10, + 'D': 1e-10, + 'E':1e-10, + 'F':1e-10, + 'G':1e-10, 'device':1e-10} results_sim = sim.run_sim('ipopt',variances=sigmas,seed=123453256, @@ -265,12 +265,12 @@ def rule_odes(m, t): dZ_Dataframe = pd.DataFrame(data=results_sim.dZdt) write_absorption_data_to_csv( - os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_S_data_input_noiselesspoints.csv'), S_Dataframe) - write_spectral_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_D_data_input_noiselesspoints.csv'), D_Dataframe) - write_concentration_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_C_data_input_noiselesspoints.csv'), C_Dataframe) - write_concentration_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_Z_data_input_noiselesspoints.csv'), Z_Dataframe) + os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_S_data_input_noiselesspoints2.csv'), S_Dataframe) + write_spectral_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_D_data_input_noiselesspoints2.csv'), D_Dataframe) + write_concentration_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_C_data_input_noiselesspoints2.csv'), C_Dataframe) + write_concentration_data_to_csv(os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_Z_data_input_noiselesspoints2.csv'), Z_Dataframe) write_concentration_data_to_csv( - os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_dZ_data_input_noiselesspoints.csv'), dZ_Dataframe) + os.path.join(dataDirectory, 'FeCaseexamplewithoutTemp_dZ_data_input_noiselesspoints2.csv'), dZ_Dataframe) if with_plots: diff --git a/kipet/examples/PaperExamples/paper_example3_inputsandnonabs.py b/kipet/examples/PaperExamples/paper_example3_inputsandnonabs.py index a91d9873..052a8333 100644 --- a/kipet/examples/PaperExamples/paper_example3_inputsandnonabs.py +++ b/kipet/examples/PaperExamples/paper_example3_inputsandnonabs.py @@ -7,7 +7,6 @@ # Sample Problem # Paper Example 3 # Inputs, extra_states and non-absorbing species -# Slight modification of the Michael's reaction # Estimation with known variances of spectral data using pyomo discretization # # C_k(t_i) = Z_k(t_i) + w(t_i) for all t_i in measurement points @@ -47,13 +46,13 @@ # components components = dict() - components['AH'] = 0.395555 + components['A'] = 0.395555 components['B'] = 0.0351202 components['C'] = 0.0 - components['BH+'] = 0.0 - components['A-'] = 0.0 - components['AC-'] = 0.0 - components['P'] = 0.0 + components['D'] = 0.0 + components['E'] = 0.0 + components['F'] = 0.0 + components['G'] = 0.0 builder.add_mixture_component(components) @@ -64,7 +63,7 @@ builder.add_algebraic_variable(algebraics) params = dict() - params['k0'] = 49.7796 + params['k0'] = 0.2545 params['k1'] = 8.93156 params['k2'] = 1.31765 params['k3'] = 0.310870 @@ -80,21 +79,21 @@ # stoichiometric coefficients gammas = dict() - gammas['AH'] = [-1, 0, 0, -1, 0] + gammas['A'] = [-1, 0, 0, -1, 0] gammas['B'] = [-1, 0, 0, 0, 1] gammas['C'] = [0, -1, 1, 0, 0] - gammas['BH+'] = [1, 0, 0, 0, -1] - gammas['A-'] = [1, -1, 1, 1, 0] - gammas['AC-'] = [0, 1, -1, -1, -1] - gammas['P'] = [0, 0, 0, 1, 1] + gammas['D'] = [1, 0, 0, 0, -1] + gammas['E'] = [1, -1, 1, 1, 0] + gammas['F'] = [0, 1, -1, -1, -1] + gammas['G'] = [0, 0, 0, 1, 1] def rule_algebraics(m, t): r = list() - r.append(m.Y[t, '0'] - m.P['k0'] * m.Z[t, 'AH'] * m.Z[t, 'B']) - r.append(m.Y[t, '1'] - m.P['k1'] * m.Z[t, 'A-'] * m.Z[t, 'C']) - r.append(m.Y[t, '2'] - m.P['k2'] * m.Z[t, 'AC-']) - r.append(m.Y[t, '3'] - m.P['k3'] * m.Z[t, 'AC-'] * m.Z[t, 'AH']) - r.append(m.Y[t, '4'] - m.P['k4'] * m.Z[t, 'AC-'] * m.Z[t, 'BH+']) + r.append(m.Y[t, '0'] - m.P['k0'] * m.Z[t, 'A'] * m.Z[t, 'B']) + r.append(m.Y[t, '1'] - m.P['k1'] * m.Z[t, 'E'] * m.Z[t, 'C']) + r.append(m.Y[t, '2'] - m.P['k2'] * m.Z[t, 'F']) + r.append(m.Y[t, '3'] - m.P['k3'] * m.Z[t, 'F'] * m.Z[t, 'A']) + r.append(m.Y[t, '4'] - m.P['k4'] * m.Z[t, 'F'] * m.Z[t, 'D']) return r #: there is no ae for Y[t,5] because step equn under rule_odes functions as the switch for the "C" equation @@ -124,12 +123,12 @@ def rule_odes(m, t): model = builder.create_pyomo_model(0, 600) - non_abs = ['AC-'] + non_abs = ['F'] builder.set_non_absorbing_species(model, non_abs) #Load data: dataDirectory = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))), 'data_sets')) - filenameD = os.path.join(dataDirectory,'FeCaseexamplewithoutTemp_D_data_input_noiselesspoints.csv') + filenameD = os.path.join(dataDirectory,'FeCaseexamplewithoutTemp_D_data_input_noiselesspoints2.csv') D_frame = read_spectral_data_from_csv(filenameD) builder.add_spectral_data(D_frame) model = builder.create_pyomo_model(0., 600.) @@ -166,10 +165,10 @@ def rule_odes(m, t): #to this function as an argument dictionary #New Inputs for discrete feeds - Z_step = {'AH': .03}#Which component and which amount is added + Z_step = {'A': .03}#Which component and which amount is added X_step = {'V': .01} jump_states = {'Z': Z_step, 'X': X_step} - jump_points1 = {'AH': 101.035}#Which component is added at which point in time + jump_points1 = {'A': 101.035}#Which component is added at which point in time jump_points2 = {'V': 303.126} jump_times = {'Z': jump_points1, 'X': jump_points2} @@ -191,19 +190,19 @@ def rule_odes(m, t): # ========================================================================= # USER INPUT SECTION - Parameter Estimation # # ======================================================================== - sigmas={'AH': 1e-10, + sigmas={'A': 1e-10, 'B':1e-10, - 'C': 1e-10, - 'BH+': 1e-10, - 'A-':1e-10, - 'P':1e-10, - 'device':1e-10} + 'C': 1e-10, + 'D': 1e-10, + 'E':1e-10, + 'G':1e-10, + 'device':1e-10} model = builder.create_pyomo_model(0., 600.)#0.51667 model.del_component(params) - builder.add_parameter('k0', init=0.9*49.7796,bounds=(0.0, 100.0)) - builder.add_parameter('k1', 8.93156) + builder.add_parameter('k0', init=0.9*0.2545,bounds=(0.0, 100.0)) + builder.add_parameter('k1', init=0.9*8.93156, bounds=(0.0,100.)) builder.add_parameter('k2', init=0.9*1.31765,bounds=(0.0,100.)) builder.add_parameter('k3', init=0.9*0.310870, bounds=(0.,100.)) builder.add_parameter('k4', init=0.9*3.87809,bounds=(0.0, 100.))