tutorial2str.py

Go to the documentation of this file.

 
 
 
import os
import sys
 
import conoptpy
import adolc
 
import numpy as np
 
sys.path.append('../../common/')
import std
 
 
class TutModelData(conoptpy.ModelData):
   def __init__(self):
      self.Al   = 0.16
      self.Ak   = 2.0
      self.Ainp = 0.16
      self.Rho  = 1.0
      self.K    = 4.0
      super().__init__()
 
   def buildModel(self):
      """
      adding the variables and constraints to the model
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      """
      # adding the variables to the model
      self.addVariable(0.1, conoptpy.CONOPT_INF, 0.5)
      self.addVariable(0.1, conoptpy.CONOPT_INF, 0.5)
      self.addVariable(0.0, conoptpy.CONOPT_INF)
      self.addVariable(0.0, conoptpy.CONOPT_INF)
 
      # adding the constraints to the model
      self.addConstraint(conoptpy.ConstraintType_Free, -0.1, [0, 1, 2, 3],
                    [-1, -1, 0, 0], [0, 0, 1, 1])
      self.addConstraint(conoptpy.ConstraintType_Eq, 0.0, [0, 1, 2], [0, 0, -1],
                    [1, 1, 0])
      self.addConstraint(conoptpy.ConstraintType_Eq, 4.0, [2, 3], [1, 2],
                    [0, 0])
 
      # setting the objective constraint
      self.setObjectiveElement(conoptpy.ObjectiveElement_Constraint, 0)
 
      # setting the optimisation direction
      self.setOptimizationSense(conoptpy.Sense_Maximize)
 
      # initialising the automatic differentiation
      self.initialiseAutoDiff()
 
      # computes and sets the Hessian structure
      self.computeHessianStructure()
 
 
   def tapeFunction(self, x, rowno):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      evaluates the nonlinear function and records a tape is necessary
 
      @param x       current point to be evaluated
      @param rowno   the index of the constraint. This is also used for the trace tag.
      """
      adolc.trace_on(rowno)
      ax = adolc.as_adouble(x)
 
      # marking the x variables as independent
      for item in iter(ax):
         item.declareIndependent()
 
      L   = ax[0]
      Inp = ax[1]
      Out = ax[2]
      P   = ax[3]
 
      ay = adolc.as_adouble(0)
      if rowno == 0:
         ay = P * Out
      elif rowno == 1:
         hold1 = (self.Al*pow(L,(-self.Rho)) + self.Ak*pow(self.K,(-self.Rho)) + self.Ainp*pow(Inp,(-self.Rho)))
         hold2 = pow(hold1,( -1./self.Rho ))
 
         ay = hold2
 
      ay.declareDependent()
      adolc.trace_off()
 
 
   def initialiseAutoDiff(self):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      initialises the automatic differentiation
      """
      x = []
      for v in range(self.numVar()):
         x.append(self.getVariable(v).curr)
 
      for c in range(self.numCons()):
         self.tapeFunction(x, c)
 
 
   def computeHessianStructure(self):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      uses the automatic differentiation methods to compute the Hessian
      structure
      """
      x = []
      for v in range(self.numVar()):
         x.append(self.getVariable(v).curr)
 
      hessstr = {}
      for c in range(self.numCons()):
         hesspat = adolc.hess_pat(c, x, 0)
 
         index = 0
         for i, n in enumerate(hesspat[0]):
            if n > 0:
               for j in range(n):
                  if i >= hesspat[1][index + j]:
                     if i not in hessstr.keys():
                        hessstr[i] = []
                     hessstr[i].append(hesspat[1][index + j])
 
            index += n
 
      rowindex = [int(r) for r in range(self.numVar()) for _ in hessstr.get(r, [])]
      colindex = [int(c) for r in range(self.numVar()) for c in hessstr.get(r, [])]
 
      # setting the structure of the second derivative of the Lagrangian
      self.setSDLagrangianStructure(rowindex, colindex)
 
 
   def evaluateNonlinearTerm(self, x, rowno, ignerr, thread):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      """
      try:
         g = adolc.function(rowno, x)[0]
      except adolc.BranchException:
         self.tapeFunction(x, rowno)
         g = adolc.function(rowno, x)[0]
 
      return g
 
 
   def evaluateNonlinearJacobian(self, x, rowno, jacnum, ignerr, thread):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      """
      jac = []
      try:
         jacres = adolc.gradient(rowno, x)
      except adolc.BranchException:
         self.tapeFunction(x, rowno)
         jacres = adolc.gradient(rowno, x)
 
      for i in jacnum:
         jac.append(jacres[i])
 
      return jac
 
 
   def evaluateSDLagrangian(self, x, u, hessianrow, hessiancol):
      """
      @ingroup PYTHON1THREAD_AD_TUTORIAL2STR
      """
      numhessian = len(hessianrow)
      hessian = [0 for i in range(numhessian)]
 
      for c in range(self.numCons()):
         try:
            hessres = adolc.hessian(c, x)
         except adolc.BranchException:
            self.tapeFunction(x, rowno)
            hessres = adolc.hessian(c, x)
 
         for i in range(numhessian):
            hessian[i] += u[c]*hessres[hessianrow[i]][hessiancol[i]]
 
      return hessian
 
if __name__ == "__main__":
   name = os.path.basename(__file__)[:-3]
 
   conopt = conoptpy.Conopt(name)
   model = TutModelData()
   msghdlr = std.TutMessageHandler(name)
 
   model.buildModel()
 
   conopt.loadModel(model)
   conopt.setMessageHandler(msghdlr)
 
   # getting the license variables
   license_int_1 = os.environ.get('CONOPT_LICENSE_INT_1', None)
   license_int_2 = os.environ.get('CONOPT_LICENSE_INT_2', None)
   license_int_3 = os.environ.get('CONOPT_LICENSE_INT_3', None)
   license_text = os.environ.get('CONOPT_LICENSE_TEXT', None)
   if license_int_1 is not None and license_int_2 is not None \
         and license_int_3 is not None and license_text is not None:
      conopt.setLicense(int(license_int_1), int(license_int_2),
                        int(license_int_3), license_text)
 
   coi_error = conopt.solve()
 
   retcode = std.checkSolve(conopt, 0.572943, coi_error)
 
   sys.exit(retcode)