DianaOptimization: OEDSigmaPoint.py

# -----------------------------------------------------------------------
#  Diana process modelling, simulation and analysis software
#  Copyright (c) 2009, Sergiy Gogolenko
#  e-mail: gogolenk@mpi-magdeburg.mpg.de
#  All rights reserved.
# -----------------------------------------------------------------------
#  Description:
#   Simple example of using DIANA for solving parameter fitting problems
#  with deterministic optimizer (in this case LBFGS-solver).

import diana, sys, os
#from math import *

modelPath = "Substrateuptake.so"
observedParams = ["c_x", "c_s"]
fileMesData = "./data/observations.dat"
integratorName = "odessa" #"ida"
optimizerPEName = "nmsimplex"
methodOEDName = "oedsigmapoint"
optimizerOEDName = "direct" #"genetic" #"nmsimplex"

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

##########################################################################
#                          PI Definition Section                         #
##########################################################################

# create optimizer for PI
sfactory = dmain.GetSolverFactory()
optimizerPE = sfactory.CreateSolver(diana.CAPE_NLP, None, optimizerPEName)

# set optimizer's parameters
pars_PEsolver = optimizerPE.GetParameters()
pars_PEsolver["StopCondition"].SetValue(1+8)
pars_PEsolver["MaxTimeOfWork"].SetValue(10)
pars_PEsolver["MaxCountOfEvals"].SetValue(1)
pars_PEsolver["VerboseLevel"].SetValue(0)

pars_PEsolver["Step"].SetValue(0.01)
pars_PEsolver["MaxFn"].SetValue(40)
pars_PEsolver["StopCr"].SetValue(1e-2)
pars_PEsolver["NLoop"].SetValue(20)
pars_PEsolver["IQuad"].SetValue(True)
pars_PEsolver["Simp"].SetValue(1e-6)

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

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

# describe estimated parameters
sps = [
    diana.DianaNLPRealParameterSpec("k_s", "", 1.0, 1.0, 3.0, 0.0001),
    diana.DianaNLPRealParameterSpec("mu_max", "",  4.0, 4.0, 6.0, 0.0001)
    ]

# create paramerter fitting task
taskPE = diana.ParameterFittingTask(dmain, integratorName, model, md, sps)
taskPE.Initialize()
print "Parameter fitting task "\
      "for %s is initialized" % model.GetComponentName()

# set prepeared PI task to PI optimizer
optimizerPE.SetNLPTask(taskPE)
optimizerPE.Initialize()
print "PI optimizer %s initialized" % optimizerPE.GetComponentName()

#########################################################################
#                          OED Definition Section                       #
#########################################################################
# describe design variables
colDesignVars = [
    diana.DianaNLPRealParameterSpec("q_0", "",  0.09, 0.07, 0.08, 0.0001)
    ]

# create OED problem
nlptaskfactory = dmain.GetNLPTaskFactory()
taskOED = nlptaskfactory.CreateOEDTask(taskPE, optimizerPE,
                                       colDesignVars, methodOEDName)
pars_OEDtask = taskOED.GetParameters()
pars_OEDtask["OptimalityCriterion"].SetValue(diana.DIANA_OED_EStar)
pars_OEDtask["ExecException"].SetValue(True)
pars_OEDtask["Alpha"].SetValue(0.5)
taskOED.Initialize()
print "OED task initialized..."

########################################################################
#                     Create OED task reporting                       #
########################################################################

# create reportOEDTask
reportOEDTask = dmain.CreateReportingInterface(methodOEDName)
# set parameters for reportOEDTask class
pars_taskOEDrep = reportOEDTask.GetParameters()
pars_taskOEDrep["OutputFolder"].SetValue("./oedtask")
pars_taskOEDrep["ConsoleOutput"].SetValue(True)
pars_taskOEDrep["SamplingPoints"].SetValue(True)

reportOEDTask.Initialize()

taskOED.SetReportingInterface(reportOEDTask)
print "OED task report initialized"

########################################################################
#                          OED Solving Section                         #
########################################################################

# create optimizer for OED problem
optimizerOED = sfactory.CreateSolver(diana.CAPE_NLP,
                                     None, optimizerOEDName)

# settings for OED optimizer
pars_OEDsolver = optimizerOED.GetParameters()
pars_OEDsolver["Eps"].SetValue(1.0000000000e-4)
pars_OEDsolver["Maxf"].SetValue(5000)
pars_OEDsolver["MaxT"].SetValue(1000)
pars_OEDsolver["Algmethod"].SetValue(False)
pars_OEDsolver["FGlPer"].SetValue(1.0000000000e-02)
pars_OEDsolver["VolPer"].SetValue(-1.0000000000e+00)
pars_OEDsolver["SigmaPer"].SetValue(-1.0000000000e+00)
pars_OEDsolver["FGlobal"].SetValue(-10.15320)
pars_OEDsolver["VerboseLevel"].SetValue(2)

# create report
report = dmain.CreateReportingInterface(optimizerOEDName)
# set parameters for report class
pars_solverOEDrep = report.GetParameters()
pars_solverOEDrep["Author"].SetValue("SG")
pars_solverOEDrep["OutputFolder"].SetValue("./oed")
pars_solverOEDrep["ConsoleOutput"].SetValue(True)
pars_solverOEDrep["Samples"].SetValue(True)

report.Initialize()
print "Report initialized"

# set prepeared task and report to OED optimizer
optimizerOED.SetReportingInterface(report)
optimizerOED.SetNLPTask(taskOED)
optimizerOED.Initialize()
print "OED optimizer %s initialized" % optimizerOED.GetComponentName()

# start optimization
print "OED by %s for %s..." % (taskOED.GetComponentName(),
                               model.GetComponentName())
try:
    optimizerOED.Solve()
except:
    print "Best solution: ", optimizerOED.GetSolution()
    raise
print "Solution: ", optimizerOED.GetSolution()
print "Parameter values:"
print task.GetSoughtParameters()