How to fit data on algebraic variables

[notesofdabbler]

For a reaction example with components A, B, C, D, the data collected at different time points is y = B/(A + B). That is a csv file with columns for time, y. In the builder, I can specify y as an algebraic variable but wasn't sure how to make the builder use the objective of sum((y_predicted - y)**2) to fit the data.

[kwmcbride]

I will have to look into this myself since I have not yet tried to solve such a problem. Do you have a basic example I can start with?

[notesofdabbler]

I have created a variation of fitting A->B->C in this notebook where measurements are on A and B/B+C. I have created the model but am not sure how to add the objective corresponding to residual of B/B+C between measured and predicted values

[kwmcbride]

Thanks for the notebook. I think the issue is going to be solved by letting the user use a custom objective when providing certain types of data, such as in your case.

[kwmcbride]

I have your example working in Kipet now. I just need to work out how to make it user friendly. A fix for this type of problem will be soon available.

[kwmcbride]

Here is the working code from the new version. This is how it will look once I make the changes to the new version, but it seems user friendly enough.

"""
Different objectives: data provided in ratios
"""
import kipet
kipet_model = kipet.KipetModel()

full_data = kipet.read_file(kipet.set_directory('ratios.csv'))

r1 = kipet_model.new_reaction('reaction-1')

# Add the model parameters
r1.add_parameter('k1', init=5.0, bounds=(0.0, 10.0))
r1.add_parameter('k2', init=5.0, bounds=(0.0, 10.0))

# Declare the components and give the initial values
r1.add_component('A', state='concentration', init=1.0)
r1.add_component('B', state='concentration', init=0.0)
r1.add_component('C', state='concentration', init=0.0)

filename = 'ratios.csv'
r1.add_dataset(data=full_data[['A']])
r1.add_dataset('y_data', category='custom', data=full_data[['y']])

# Define the reaction model
def rule_odes(m,t):
    exprs = dict()
    exprs['A'] = -m.P['k1']*m.Z[t,'A']
    exprs['B'] = m.P['k1']*m.Z[t,'A']-m.P['k2']*m.Z[t,'B']
    exprs['C'] = m.P['k2']*m.Z[t,'B']
    return exprs 

r1.add_equations(rule_odes)

# This will let you plot the results against the given ratio B/(B + C)
r1.add_algebraic_variables('y', init = 0.0, bounds = (0.0, 1.0))

def rule_algebraics(m, t):
    r = list()
    r.append(m.Y[t, 'y']*(m.Z[t, 'B'] + m.Z[t, 'C']) - m.Z[t, 'B'])
    return r

r1.add_algebraics(rule_algebraics)

def rule_objective(m, t):
    return ((m.UD[t, 'y']*(m.Z[t, 'B'] + m.Z[t, 'C'] + 1e-12) - m.Z[t, 'B'])**2)

r1.add_custom_objective(rule_objective)   

solve = True

if solve:

    r1.run_opt()
    # Display the results
    r1.results.show_parameters

    # New plotting methods
    r1.results.plot()

And the results (please ignore the labelling on the second image!):

[notesofdabbler]

This is great and looks quite user friendly. Would this approach also work in the case where the time points where data for 'A' is available is not the same as the time points where data for 'y' is available (like the heterogeneous data example)

[kwmcbride]

I made a small change to the TemplateBuilder to make sure that this works. So, yes, you will be able to use data with mixed time points.

I would recommend that you use enough concentration points along with the ratios. If you only use ratios, you won't be able to solve it.

[kwmcbride]

I changed it to where you only need to declare the algebraic rule and then state that you want to use an algebraic variable as an objective. It needs to match the custom data that you enter. It will match the data with the variable using least squares terms added to the objective function.

r1.add_custom_objective('y')

See the following complete example:

import kipet

kipet_model = kipet.KipetModel()

full_data = kipet.read_file(kipet.set_directory('ratios.csv'))

r1 = kipet_model.new_reaction('reaction-1')

# Add the model parameters
r1.add_parameter('k1', init=5.0, bounds=(0.0, 10.0))
r1.add_parameter('k2', init=5.0, bounds=(0.0, 10.0))

# Declare the components and give the initial values
r1.add_component('A', state='concentration', init=1.0)
r1.add_component('B', state='concentration', init=0.0)
r1.add_component('C', state='concentration', init=0.0)

filename = 'ratios.csv'

r1.add_dataset(data=full_data[['A']], remove_negatives=True)
r1.add_dataset('y_data', category='custom', data=full_data[['y']])

# Define the reaction model
def rule_odes(m,t):
    exprs = dict()
    exprs['A'] = -m.P['k1']*m.Z[t,'A']
    exprs['B'] = m.P['k1']*m.Z[t,'A']-m.P['k2']*m.Z[t,'B']
    exprs['C'] = m.P['k2']*m.Z[t,'B']
    return exprs 

r1.add_equations(rule_odes)
r1.add_algebraic_variables('y', init = 0.0, bounds = (0.0, 1.0))

def rule_algebraics(m, t):
    r = list()
    r.append(m.Y[t, 'y']*(m.Z[t, 'B'] + m.Z[t, 'C'] + 1e-20) - m.Z[t, 'B'])
    return r

r1.add_algebraics(rule_algebraics)
r1.add_custom_objective('y')

r1.run_opt()
r1.results.show_parameters
r1.results.plot()

[notesofdabbler]

thanks. Will this be part of the new pull request called "new data" that I saw in the repo?

[kwmcbride]

At the moment this is not part of new_data. I have two further updates to KIPET after new_data, one containing fixes for Pyomo 5.7 and the update above is based on that version. Once new_data is approved (I have to finish the documentation), this change should be quickly accepted.

[kwmcbride]

I have added this to the pull request "new_data".