How to get this Parameter Estimation code to run?

[Manuel]

Dear Dragan,

again, I am trying to do optimization with DAE Tools and find it a bit difficult to get it to work.
I am able to run the model with input values taken from a CSV file.

python testODM2prod_ParEst.py simulation

Now, I have measurement values (also in a CSV file), and I would like to adjust the model's parameters, so the simulated states fit to the measured states.
First, I tried using the steps from opt_tutorial5, but that didn't work. It seems that your independent variable is x, why mine is t. Therefore, I cannot use ReAssignValue in the loop in Function.
Second, I tried using the steps from tutorial_che_opt_4; I got a version that would execute ( via python testODM2prod_ParEst.py console, see below), but the results need to be vastly improved. One reason is that I use a custom Run() method for simulation, and not sure how to merge it with the Run() method for optimization. Another reason is that one of the output variables I am trying to fit is algebraic, so the assumption that L2(t=0) = 0 might be False.
Maybe you could take a look at the code and tell me how to get it to work.

Kind regards,
Manuel

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
ODM2prod in Python!!!
@author: rotton
"""

__doc__ = """
Simulating ODM2prod in DAE Tools; this time with real feeding
"""

import numpy, os, sys, tempfile
from time import localtime, strftime
from daetools.solvers.ipopt import pyIPOPT
from daetools.pyDAE import *  #bad style, but who cares?
from daetools.pyDAE.data_reporters import daeCSVFileDataReporter

from pyUnits import m, kg, s, K, bar, mol, l, kmol

from daetools.pyDAE.variable_types import volume_t, mass_concentration_t, volume_flowrate_t
#Create custom variable types
rate_constant_t =  daeVariableType("rate_constant_t", s**(-1), 1e-8, 1e3, 6e-6, 1e-7)
fummel_t =  daeVariableType("fummel_t", (m**3/s)*(m**3/kg), 1e-8, 1e4, 7e-4, 1e-6)
mass_reaction_rate_t =  daeVariableType("mass_reaction_rate_t", kg/m**3/s, -100, 1e10, 1e-5, 1e-7)
L2_t = daeVariableType("L2_t", unit(), -1.0e+20, 1.0e+20, 0.0, 1e-07)

class anaerobicDigestionModel(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)
        # Model parameters --> as variables for PI/optimization!
        self.k_1 = daeVariable("k_1", rate_constant_t, self, "ODM1 decay rate")
        self.f_1 = daeVariable("f_1", fummel_t, self, "Sub2prod conversion factor")
        # Physical parameters
        self.V_liq = daeParameter("V_{liq}", m**3, self, "Volume of the liquid phase")
        # Differential variables of model
        self.X = daeVariable("ODM1", mass_concentration_t, self, \
                             "ODM1 concentration (in the liquid phase)")
        self.Y = daeVariable("ODM2", mass_concentration_t, self, \
                             "ODM2 concentration (in the liquid phase)")
        # Inlet variable(s)
        self.q_in = daeVariable("q_{in}", volume_flowrate_t, self, "Volumetric flowrate_inlet")
        self.X_in = daeVariable("ODM1_{in}", mass_concentration_t, self, \
                                    "ODM1 concentration in the liquid phase_inlet")
        self.Y_in = daeVariable("ODM2_{in}", mass_concentration_t, self, \
                                    "ODM2 concentration in the liquid phase_inlet")
        # Algebraic variable
        self.q_prod = daeVariable("q_{prod}", volume_flowrate_t, self, \
                                  "Volumetric flowrate of product")

        # Reaction rate, defined as separate algebraic variable
        self.rho_1 = daeVariable("ρ_1", mass_reaction_rate_t, self, "ODM1 reaction rate")

    def DeclareEquations(self):
        # Create adouble objects to make equations more readable
        X = self.X()
        Y = self.Y()
        q_prod = self.q_prod()
        k_1 = self.k_1()
        f_1 = self.f_1()
        # Reaction rates
        rho_1 = self.rho_1()
        # Inlet variables
        q_in = self.q_in()
        X_in = self.X_in()
        Y_in = self.Y_in()
        # Physical parameters
        V_liq = self.V_liq()

        daeModel.DeclareEquations(self)
        #Differential equations
        eq = self.CreateEquation("MassBalance01", "Mass balance for ODM1")
#        eq.BuildJacobianExpressions = True
        eq.Residual = q_in/V_liq*(X_in - X) - rho_1 - dt(X)

        eq = self.CreateEquation("MassBalance02", "Mass balance for ODM2")
        eq.Residual = q_in/V_liq*(Y_in - Y) - dt(Y)

        #Algebraic equations
        eq = self.CreateEquation("MassBalance03", "Product formation")
        eq.Residual = 0.52*f_1*X - q_prod

        # Reaction rates
        eq = self.CreateEquation("Rate01", "ODM reaction rate")
        eq.Residual = k_1*X - rho_1

class anaerobicDigestionSimulation(daeSimulation):
    def __init__(self):
        daeSimulation.__init__(self)

        self.m = anaerobicDigestionModel("ODM_2h")
        self.m.Description = __doc__

    def SetUpParametersAndDomains(self):
        # Physical parameters
        self.m.V_liq.SetValue(159 * m**3)

    def SetUpVariables(self):
        #X_steady = q_in/V_liq*X_in/(q_in/V_liq + k_1) = X_0
        self.m.X.SetInitialCondition(11.55 * kg/m**3)
        self.m.Y.SetInitialCondition(11.55*0.2 * kg/m**3)
        #Initial value for variable receiving external values
        self.m.q_in.AssignValue(0.0)
        self.m.X_in.AssignValue(0.0)
        self.m.Y_in.AssignValue(0.0)
        #Parameters --> as variables for PI/optimization!
        self.m.k_1.AssignValue(0.19/86400)
        self.m.f_1.AssignValue(29.0/86400)

    def Run(self):
        self.integratorStats = {}

        # Create a reversed list to serve as a lifo stack.
        feedingSchedule_stack = list(feedingSchedule)
        feedingSchedule_stack.reverse()

        # Check the last feed-in time (it has to be less than the specified TimeHorizon)
        if feedingSchedule_stack:
            lastFeedingTime = feedingSchedule_stack[0][0]
            if lastFeedingTime > self.TimeHorizon:
                raise RuntimeWarning('The last feed-in time is past the specified time horizon')

        while self.CurrentTime < self.TimeHorizon:
            # Get the next reporting time (using self.NextReportingTime which is based on
            # TimeHorizon and ReportingInterval). Do not allow to get past the TimeHorizon
            # (this should be handled by daetools, but double-checking is fine).
            t_rep = self.NextReportingTime
            if t_rep > self.TimeHorizon:
                t_rep = self.TimeHorizon

            self.Log.Message('Integrating from %f to %f ...' % (self.CurrentTime, t_rep), 0)
            while True:
                # If the feeding schedule is empty (or it has already been consumed) break the while loop.
                if len(feedingSchedule_stack) == 0:
                    break

                # Get the feeding data from the stack (located in its tail - the last item)
                t_feed, q_in, X_in, Y_in = feedingSchedule_stack[-1]
                # If the next feeding moment belongs to this reporting interval it has to be 
                # processed. Otherwise, break this loop and simply integrate to the next
                # reporting time.
                if t_feed > t_rep:
                    break

                # The next feed-in moment belongs to this time interval - remove it from the stack.
                feedingSchedule_stack.pop()

                # Integrate to the specified feed-in time, report data and set progress.
                self.Log.Message('  Integrating to the next feed-in time %8.1f...' % (t_feed), 0)
                # set reportDataAroundDiscontinuities argument to False to avoid reporting around
                # discontinuities. Anyway, you do not have discontinuities in your model so it
                # doesn't really matter.
                self.IntegrateUntilTime(t_feed, eDoNotStopAtDiscontinuity, False) 
                # Report the results only for debugging purposes (to overlook what happens
                # during simulation).

                self.integratorStats[self.CurrentTime] = self.DAESolver.IntegratorStats

                # Once we reached the next feed-in time, re-set q_in and A_in and re-initialise
                # the system.
                self.m.q_in.ReAssignValue(q_in/86400) # values in CSV file are in m**3/day
                self.m.X_in.ReAssignValue(X_in)
                self.m.Y_in.ReAssignValue(Y_in)
                self.Reinitialize()

            # Finalise the current interval:
            #   If the target reporting time has not been reached, integrate to t_rep, report the
            # data and set progress.
            if self.CurrentTime < t_rep:
                self.Log.Message('  Finishing the current reporting interval (integrate to %f)...'\
                                 % (t_rep), 0)
                self.IntegrateUntilTime(t_rep, eDoNotStopAtDiscontinuity, False)
                self.ReportData(self.CurrentTime)
                self.Log.SetProgress(int(100.0 * self.CurrentTime/self.TimeHorizon))


        self.Log.Message('Simulation has finished successfully (the final time reached: %f)' % \
                         self.CurrentTime, 0)
        totalDAESteps = 0
        for t, stats in self.integratorStats.items():
            totalDAESteps += stats['NumSteps']

        print('\n\nIntegrators stats:')
        print('Total number of DAE steps = %d' % totalDAESteps)

#   anaerobicDigestionModel_Opt inherits all parameters/variables from the base class.
class anaerobicDigestionModel_Opt(anaerobicDigestionModel):
    def __init__(self, Name, Parent = None, Description = ""):
        anaerobicDigestionModel.__init__(self, Name, Parent, Description)

        # Observed values at the specific time interval
        self.X_obs = daeVariable("X_obs", no_t, self, \
                                 "Observed value 1 at the specified time interval")
        self.q_obs = daeVariable("q_obs", no_t, self, \
                                 "Observed value 2 at the specified time interval")

        # This L2 norm sums all L2 norms in the previous time intervals
        self.L2      = daeVariable("L2",      L2_t, self, \
                                   "Current L2 norm: ||yi(t) - yi_obs(t)||^2")
        self.L2_prev = daeVariable("L2_prev", L2_t, self, "L2 norm in previous time intrvals")

    def DeclareEquations(self):
        anaerobicDigestionModel.DeclareEquations(self)

        # L2-norm ||yi(t) - yi_obs(t)||^2
        # L2 norm is a sum of the L2 norm in the previous time steps (L2_prev)
        # and the current norm: s1 + s2.
        # L2_prev will be reset after every time interval where we have observed values.
        s1 = (self.X() - self.X_obs())**2
        s2 = (self.q_prod() - self.q_obs())**2
        eq = self.CreateEquation("L2", "")
        eq.Residual = self.L2() - (self.L2_prev() + s1 + s2)
        eq.CheckUnitsConsistency = False

# Simulation class that will be used by the optimisation.
class anaerobicDigestionSimulation_opt(daeSimulation):
    def __init__(self):
        daeSimulation.__init__(self)
        self.m = anaerobicDigestionModel_Opt("Odm2prod_ParEst")
        self.m.Description = __doc__

    def SetUpParametersAndDomains(self):
        # Physical parameters
        self.m.V_liq.SetValue(159 * m**3)

    def SetUpVariables(self):
        #Parameters --> as variables for PI/optimization!
        self.m.k_1.AssignValue(0.19/86400)
        self.m.f_1.AssignValue(29.0/86400)

        self.m.X.SetInitialCondition(X_t0 * kg/m**3)
        self.m.Y.SetInitialCondition(X_t0*0.2 * kg/m**3)
        #Initial value for variable receiving external values
        self.m.q_in.AssignValue(0.0)
        self.m.X_in.AssignValue(0.0)
        self.m.Y_in.AssignValue(0.0)

        # Initialise variables required for parameter estimation.
        # Notate bene: Observed values should match initial conditions at t = 0
        # --> Not the case here, since q_prod is algebraic and depends directly on X
        self.m.X_obs.AssignValue(X_t0)
        self.m.q_obs.AssignValue(q_obs_t0)
        self.m.L2_prev.AssignValue(0.0)

    def Run(self):
        for t, tn in enumerate(times):
            # Reset L2_prev value to the current L2
            if t == 0:
                self.m.L2_prev.ReAssignValue(0.0)
            else:
                L2 = self.m.L2.GetValue()
                self.m.L2_prev.ReAssignValue(L2)

            # Reset observed values to match the current interval end time
            self.m.X_obs.ReAssignValue(t_xq_obs[t+1,1])
            self.m.q_obs.ReAssignValue(t_xq_obs[t+1,2]/3600*0.52) #q_obs in CSV file in m^3/h!

            # Reinitialise the DAE system after all changes made above
            self.Reinitialize()

            # Integrate, report data and set progress
            self.Log.Message('Integrating from %f to %f ...' % (self.CurrentTime, tn), 0)
            self.IntegrateUntilTime(tn, eDoNotStopAtDiscontinuity)
            self.ReportData(self.CurrentTime)
            self.Log.SetProgress(int(100.0 * self.CurrentTime/self.TimeHorizon))

    def SetUpOptimization(self):
        # Minimise L2-norm ||yi(t) - yi_obs(t)||^2
        self.ObjectiveFunction.Residual = self.m.L2()

        self.SetContinuousOptimizationVariable(self.m.k_1, 1e-9, 10,  0.19/86400)
        self.SetContinuousOptimizationVariable(self.m.f_1, 1e-9, 1e3, 29.0/86400)

# Read input schedule and measurement data from CSV file
feedingSchedule = numpy.genfromtxt('Inputarray_Odm2prod_2h_7d.csv', delimiter=',')
t_xq_obs = numpy.genfromtxt('Measurements_1h_7d.csv', delimiter=',')
times       = t_xq_obs[1:,0] #exclude 0
# Initial conditions
X_t0        = t_xq_obs[0, 1]
q_obs_t0    = t_xq_obs[0, 2]/3600*0.52  #q_obs in CSV file in m^3/h!

def setOptions(nlpsolver):
    nlpsolver.SetOption('print_level', 5)
    nlpsolver.SetOption('tol', 1e-6)
    #nlpsolver.SetOption('mu_strategy', 'adaptive')
    #nlpsolver.SetOption('obj_scaling_factor', 10.0)
    nlpsolver.SetOption('nlp_scaling_method', 'none') #'user-scaling')

# Setup everything manually and run in a console
def consoleSimulation():
    # First change the config values, overwriting Core/Activity/IDAS parameters
    cfg = daeGetConfig()
#    cfg.SetBoolean("daetools.core.printInfo", True)
#    cfg.SetBoolean("daetools.IDAS.printInfo", True)

    # Then create objects
    log          = daePythonStdOutLog()
    daesolver    = daeIDAS()
    datareporter = daeCSVFileDataReporter() #daeTCPIPDataReporter()
    simulation   = anaerobicDigestionSimulation()

    simulation.TimeHorizon       = 7*86400
    simulation.ReportingInterval = 1*3600
    #Specify variables which should be written to output file (all or only selected ones)
    simulation.m.SetReportingOn(True)
    simulation.ReportTimeDerivatives = False

    # Connect data reporter
    modelName = simulation.m.Name
    simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())
    directory = tempfile.gettempdir()
    csv_filename  = os.path.join(directory, "%s.csv"  % modelName)
    if(datareporter.Connect(csv_filename, simName) == False):
        sys.exit()

    # Initialize the simulation
    simulation.Initialize(daesolver, datareporter, log)

    # Save the model report and the runtime model report
    simulation.m.SaveModelReport(simulation.m.Name + ".xml")
    simulation.m.SaveRuntimeModelReport(simulation.m.Name + "-rt.xml")

    #Solve at time = 0 (initialisation)
    simulation.SolveInitial()

    # Run
    simulation.Run()
    simulation.Finalize()

def run(**kwargs):
    simulation = anaerobicDigestionSimulation_opt()
    nlpsolver  = pyIPOPT.daeIPOPT()
    relativeTolerance = 1e-6
    reportingTimes = times.tolist()
    return daeActivity.optimize(simulation, reportingInterval       = 1*3600, 
                                            timeHorizon             = 7*86400,
                                            reportingTimes          = reportingTimes,
                                            nlpsolver               = nlpsolver,
                                            nlpsolver_setoptions_fn = setOptions,
                                            relativeTolerance       = relativeTolerance,
                                            **kwargs)

if __name__ == "__main__":
    if len(sys.argv) > 1 and (sys.argv[1] == 'simulation'):
        consoleSimulation()
    else:
        guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True
        run(guiRun = guiRun)

[ciroki]

Hi Manuel,

You are right, there are problems with the optimisation part. In addition to a number of people already reported it, I am aware of that too. The optimisation support has not seen an update for many years so a major revision is required. This is true in particuar for IPOPT which has a new interface.

I'll have a look and let you know what can be done.

Dragan

[Manuel]

Dear Dragan,
thanks for looking into this.
Actually, I think my optimization problem could be implemented with the current tools. I just don't know how. My biggest question mark is how to combine both Run() methods. Maybe you could give me a hint on that. On one hand, I need the discontinuity handling, on the other hand, I need to properly address the objective function calculation.
Regards,
Manuel

[Manuel]

I've assembled an easier version, with q_in = 0. Also, I am ignoring units, so I can use days instead of seconds as base time unit.
Unfortunately, integration does no start, because of failed initilization:

Sundials IDAS solver cowardly failed to calculate initial conditions at TIME = 0; IDA_NO_RECOVERY
IDAS solver error in module 'IDAS' in function 'IDACalcIC': The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover. [IDA_NO_RECOVERY]

Even setting daetools.core.printInfo and daetools.IDAS.printInfo to True does not display additional information.
Could you provide help to find the error?

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

__doc__ = """
Simulating ODM2prod in DAE Tools; this time without feeding
"""

import numpy, os, sys, tempfile
from time import localtime, strftime
from daetools.solvers.ipopt import pyIPOPT
from daetools.pyDAE import *  #bad style, but who cares?
from daetools.pyDAE.data_reporters import daeCSVFileDataReporter

from pyUnits import m, kg, s

from daetools.pyDAE.variable_types import no_t
L2_t = daeVariableType("L2_t", unit(), -1.0e+20, 1.0e+20, 0.0, 1e-07)

class anaerobicDigestionModel(daeModel):
    def __init__(self, Name, Parent = None, Description = ""):
        daeModel.__init__(self, Name, Parent, Description)
        # Model parameters --> as variables for PI/optimization!
        self.k_1 = daeVariable("k_1", no_t, self, "ODM1 decay rate")
        self.f_1 = daeVariable("f_1", no_t, self, "Sub2prod conversion factor")
        # Physical parameters
        self.V_liq = daeParameter("V_{liq}", dimless, self, "Volume of the liquid phase")
        # Differential variables of model
        self.X = daeVariable("ODM1", no_t, self, \
                             "ODM1 concentration (in the liquid phase)")
        self.Y = daeVariable("ODM2", no_t, self, \
                             "ODM2 concentration (in the liquid phase)")
        # Inlet variable(s)
        self.q_in = daeVariable("q_{in}", no_t, self, "Volumetric flowrate_inlet")
        self.X_in = daeVariable("ODM1_{in}", no_t, self, \
                                    "ODM1 concentration in the liquid phase_inlet")
        self.Y_in = daeVariable("ODM2_{in}", no_t, self, \
                                    "ODM2 concentration in the liquid phase_inlet")
        # Algebraic variable
        self.q_prod = daeVariable("q_{prod}", no_t, self, \
                                  "Volumetric flowrate of product")

        # Reaction rate, defined as separate algebraic variable
        self.rho_1 = daeVariable("ρ_1", no_t, self, "ODM1 reaction rate")

    def DeclareEquations(self):
        # Create adouble objects to make equations more readable
        X = self.X()
        Y = self.Y()
        q_prod = self.q_prod()
        k_1 = self.k_1()
        f_1 = self.f_1()
        # Reaction rates
        rho_1 = self.rho_1()
        # Inlet variables
        q_in = self.q_in()
        X_in = self.X_in()
        Y_in = self.Y_in()
        # Physical parameters
        V_liq = self.V_liq()

        daeModel.DeclareEquations(self)
        #Differential equations
        eq = self.CreateEquation("MassBalance01", "Mass balance for ODM1")
#        eq.BuildJacobianExpressions = True
        eq.Residual = q_in/V_liq*(X_in - X) - rho_1 - dt(X)
        eq.CheckUnitsConsistency = False

        eq = self.CreateEquation("MassBalance02", "Mass balance for ODM2")
        eq.Residual = q_in/V_liq*(Y_in - Y) - dt(Y)
        eq.CheckUnitsConsistency = False

        #Algebraic equations
        eq = self.CreateEquation("MassBalance03", "Product formation")
        eq.Residual = 0.52*f_1*X - q_prod

        # Reaction rates
        eq = self.CreateEquation("Rate01", "ODM reaction rate")
        eq.Residual = k_1*X - rho_1

class anaerobicDigestionSimulation(daeSimulation):
    def __init__(self):
        daeSimulation.__init__(self)

        self.m = anaerobicDigestionModel("ODM_Batch_18h")
        self.m.Description = __doc__

    def SetUpParametersAndDomains(self):
        # Physical parameters
        self.m.V_liq.SetValue(159)

    def SetUpVariables(self):
        #X_steady = q_in/V_liq*X_in/(q_in/V_liq + k_1) = X_0
        self.m.X.SetInitialCondition(80.0)
        self.m.Y.SetInitialCondition(16.0)
        #Initial value for variable receiving external values
        self.m.q_in.AssignValue(0.0)
        self.m.X_in.AssignValue(0.0)
        self.m.Y_in.AssignValue(0.0)
        #Parameters --> as variables for PI/optimization!
        self.m.k_1.AssignValue(0.8296)
        self.m.f_1.AssignValue(9.3589)

#   anaerobicDigestionModel_Opt inherits all parameters/variables from the base class.
class anaerobicDigestionModel_Opt(anaerobicDigestionModel):
    def __init__(self, Name, Parent = None, Description = ""):
        anaerobicDigestionModel.__init__(self, Name, Parent, Description)

        # Observed values at the specific time interval
#        self.X_obs = daeVariable("X_obs", no_t, self, \
#                                 "Observed value 1 at the specified time interval")
        self.q_obs = daeVariable("q_obs", no_t, self, \
                                 "Observed value 2 at the specified time interval")

        # This L2 norm sums all L2 norms in the previous time intervals
        self.L2      = daeVariable("L2",      L2_t, self, \
                                   "Current L2 norm: ||yi(t) - yi_obs(t)||^2")
        self.L2_prev = daeVariable("L2_prev", L2_t, self, "L2 norm in previous time intrvals")

    def DeclareEquations(self):
        anaerobicDigestionModel.DeclareEquations(self)

        # L2-norm ||yi(t) - yi_obs(t)||^2
        # L2 norm is a sum of the L2 norm in the previous time steps (L2_prev)
        # and the current norm: s1 + s2.
        # L2_prev will be reset after every time interval where we have observed values.
#        s1 = (self.X() - self.X_obs())**2
        s2 = (self.q_prod() - self.q_obs())**2
        eq = self.CreateEquation("L2", "")
        eq.Residual = self.L2() - (self.L2_prev() + s2)
        eq.CheckUnitsConsistency = False

# Simulation class that will be used by the optimisation.
class anaerobicDigestionSimulation_opt(daeSimulation):
    def __init__(self):
        daeSimulation.__init__(self)
        self.m = anaerobicDigestionModel_Opt("Odm2prod_BatchParEst")
        self.m.Description = __doc__

    def SetUpParametersAndDomains(self):
        # Physical parameters
        self.m.V_liq.SetValue(159)

    def SetUpVariables(self):
        #Parameters --> as variables for PI/optimization!
        self.m.k_1.AssignValue(0.8296)
        self.m.f_1.AssignValue(9.3589)

        self.m.X.SetInitialCondition(80.0)
        self.m.Y.SetInitialCondition(16.0)
        #Initial value for variable receiving external values
        self.m.q_in.AssignValue(0.0)
        self.m.X_in.AssignValue(0.0)
        self.m.Y_in.AssignValue(0.0)

        # Initialise variables required for parameter estimation.
        # Notate bene: Observed values should match initial conditions at t = 0
        # --> Not the case here, since q_prod is algebraic and depends directly on X
#        self.m.X_obs.AssignValue(X_t0)
        self.m.q_obs.AssignValue(q_obs_t0)
        self.m.L2_prev.AssignValue(0.0)

    def Run(self):
        for t, tn in enumerate(times):
            # Reset L2_prev value to the current L2
            if t == 0:
                self.m.L2_prev.ReAssignValue(0.0)
            else:
                L2 = self.m.L2.GetValue()
                self.m.L2_prev.ReAssignValue(L2)

            # Reset observed values to match the current interval end time
#            self.m.y1_obs.ReAssignValue(y1_obs[t])
            self.m.q_obs.ReAssignValue(q_obs[t])

            # Reinitialise the DAE system after all changes made above
            self.Reinitialize()

            # Integrate, report data and set progress
            self.Log.Message('Integrating from %f to %f ...' % (self.CurrentTime, tn), 0)
            self.IntegrateUntilTime(tn, eDoNotStopAtDiscontinuity)
            self.ReportData(self.CurrentTime)
            self.Log.SetProgress(int(100.0 * self.CurrentTime/self.TimeHorizon))

    def SetUpOptimization(self):
        # Minimise L2-norm ||yi(t) - yi_obs(t)||^2
        self.ObjectiveFunction.Residual = self.m.L2()

        self.SetContinuousOptimizationVariable(self.m.k_1, 1e-9, 10,  0.8296)
        self.SetContinuousOptimizationVariable(self.m.f_1, 1e-9, 1e3, 9.3589)

# Read input schedule and measurement data from CSV file
t_xq_obs = numpy.genfromtxt('MeasuredAlgebrProductionRate_30min_18h.csv', delimiter=';')
# Initial conditions
q_obs_t0    = 20.0*24
times       = t_xq_obs[:, 0];           # times in CSV file in days!
q_obs       = t_xq_obs[:, 2]*24*0.52    # q_obs in CSV file in m^3/h!

def setOptions(nlpsolver):
    nlpsolver.SetOption('print_level', 5)
    nlpsolver.SetOption('tol', 1e-6)
    #nlpsolver.SetOption('mu_strategy', 'adaptive')
    #nlpsolver.SetOption('obj_scaling_factor', 10.0)
    nlpsolver.SetOption('nlp_scaling_method', 'none') #'user-scaling')

# Setup everything manually and run in a console
def consoleSimulation():
    # First change the config values, overwriting Core/Activity/IDAS parameters
    cfg = daeGetConfig()
    cfg.SetBoolean("daetools.core.printInfo", True)
    cfg.SetBoolean("daetools.IDAS.printInfo", True)
    cfg.SetBoolean("daetools.IDAS.LineSearchOffIC", False)

    # Then create objects
    log          = daePythonStdOutLog()
    daesolver    = daeIDAS()
    datareporter = daeCSVFileDataReporter() #daeTCPIPDataReporter()
    simulation   = anaerobicDigestionSimulation()

    simulation.TimeHorizon       = 18/24
    simulation.ReportingInterval = (1/2)/24
    #Specify variables which should be written to output file (all or only selected ones)
    simulation.m.SetReportingOn(True)
    simulation.ReportTimeDerivatives = False

    # Connect data reporter
    modelName = simulation.m.Name
    simName = simulation.m.Name + strftime(" [%d.%m.%Y %H:%M:%S]", localtime())
    directory = tempfile.gettempdir()
    csv_filename  = os.path.join(directory, "%s.csv"  % modelName)
    if(datareporter.Connect(csv_filename, simName) == False):
        sys.exit()

    # Initialize the simulation
    simulation.Initialize(daesolver, datareporter, log)

    # Save the model report and the runtime model report
    simulation.m.SaveModelReport(simulation.m.Name + ".xml")
    simulation.m.SaveRuntimeModelReport(simulation.m.Name + "-rt.xml")

    #Solve at time = 0 (initialisation)
    simulation.SolveInitial()

    # Run
    simulation.Run()
    simulation.Finalize()

def run(**kwargs):
    simulation = anaerobicDigestionSimulation_opt()
    nlpsolver  = pyIPOPT.daeIPOPT()
    relativeTolerance = 1e-6
    reportingTimes = times.tolist()
    return daeActivity.optimize(simulation, reportingInterval       = (1/2)/24,
                                            timeHorizon             = 18/24,
                                            reportingTimes          = reportingTimes,
                                            nlpsolver               = nlpsolver,
                                            nlpsolver_setoptions_fn = setOptions,
                                            relativeTolerance       = relativeTolerance,
                                            **kwargs)

if __name__ == "__main__":
    if len(sys.argv) > 1 and (sys.argv[1] == 'simulation'):
        consoleSimulation()
    else:
        guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True
        run(guiRun = guiRun)