Skip to content

Commit

Permalink
Merge pull request #100 from salvadorgarciamunoz/paperexample3new
Browse files Browse the repository at this point in the history 
Updated example 3
  • Loading branch information
@mchlshort
mchlshort authored Jun 3, 2019
2 parents a59da47 + a7d9785 commit 3a703c6
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 60 deletions.
62 changes: 31 additions & 31 deletions kipet/examples/PaperExamples/paper_example3_generatesimabsdata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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.)

Expand Down Expand Up @@ -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}

Expand All @@ -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,
Expand All @@ -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:

Expand Down
57 changes: 28 additions & 29 deletions kipet/examples/PaperExamples/paper_example3_inputsandnonabs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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.)
Expand Down Expand Up @@ -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}

Expand All @@ -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.))
Expand Down

0 comments on commit 3a703c6

Please sign in to comment.