DianaOptimization: SensParEstimation.py

# ---------------------------------------------------------------------
#  Diana process modelling, simulation and analysis software
#  Copyright (c) 2006, Sergiy Gogolenko
#  e-mail: gogolenk@mpi-magdeburg.mpg.de
#  All rights reserved.
# ---------------------------------------------------------------------
#  $Id: HafkeReactor.py Fri Aug 11 13:22:08 2006 gogolenk $
# ---------------------------------------------------------------------
#  Written under: i586-suse-linux 
# ---------------------------------------------------------------------
#  Last update: Fri Aug 11 13:22:08 2006 
# ---------------------------------------------------------------------
#  Description:
#   Simple example of using DIANA for solving parmeter fitting problems
#  with deterministic optimizer (in this case LBFGS-solver).

import diana, sys, os

modelPath = "Substrateuptake.so"
observedParams = ["c_x", "c_s"]
fileMesData = "./data/observations.dat"
integratorName = "ida" # "odessa"
optimizerName = "lbfgsb" # "ipopt"

# initialize Diana main class
dmain=diana.GetDianaMain(sys.argv);

# create NMSimplex-optimizer
sfactory = dmain.GetSolverFactory();
solver = sfactory.CreateSolver(diana.CAPE_NLP, None, optimizerName);

# set solver's parameters L-BFGS-B
pars_PEsolver = solver.GetParameters()
pars_PEsolver["StopCondition"].SetValue(1+8)
pars_PEsolver["MaxTimeOfWork"].SetValue(200)
pars_PEsolver["MaxCountOfEvals"].SetValue(1000)
pars_PEsolver["MaxCorrections"].SetValue(100)
pars_PEsolver["VerboseLevel"].SetValue(0)
pars_PEsolver["Factr"].SetValue(1e+7)
pars_PEsolver["PGradTol"].SetValue(1e-5)

# create model
mmanager=dmain.GetModelManager()
model = mmanager.CreateModel(diana.CAPE_CONTINUOUS, modelPath)
model.Initialize()

# load measured data (state values and time points) from specified file
md = diana.DianaMeasuredData(model, observedParams)
md.load(fileMesData)

# create paramerter fitting task
taskPE = diana.SensParameterFittingTask(dmain, integratorName, model, md)

# describe estimated parameters
taskPE.AddEstimatedParameterBySpec(
    diana.DianaNLPRealParameterSpec("mu_max", "", 4.0, 4.0, 6.0, 0.0001),
    diana.DIANA_PARAMETER)
taskPE.AddEstimatedParameterBySpec(
    diana.DianaNLPRealParameterSpec("k_s", "", 3.0, 1.0, 3.0, 0.0001),
    diana.DIANA_PARAMETER) 

taskPE.Initialize()
print "Parameter fitting task for %s is initialized" \
      % model.GetComponentName()

# create report
report = dmain.CreateReportingInterface(optimizerName)
# Set parameters for report class
pars_rep = report.GetParameters()
pars_rep["Author"].SetValue("SG")
pars_rep["OutputFolder"].SetValue("./parfit")
pars_rep["ConsoleOutput"].SetValue(True)
pars_rep["Samples"].SetValue(True)

report.Initialize()
print "Report initialized"

# set prepeared task and report to solver
solver.SetReportingInterface(report);
solver.SetNLPTask(taskPE);
solver.Initialize();
print "Optimizer %s initialized" % solver.GetComponentName()

# start optimization
print "Parameter fitting for %s..." % model.GetComponentName()
try:
    solver.Solve()
except:
    print "Best solution: ", solver.GetSolution()
    raise

print "Solution: ", solver.GetSolution()
print "Parameter values:"
print taskPE.GetSoughtParameters()