qp3.py

Go to the documentation of this file.

 
import os
import sys
 
import pyconopt
 
sys.path.append('../common/')
import std
 
 
class QPModelData(pyconopt.ModelData):
   def __init__(self):
      self.NN = 1000
      self.NQ = self.NN*2 - 1
      # only the lower triangle of the Q matrix is stored, since it is
      # symmetric. The diagonal has the value 10. While the first off diagonal
      # has the value 0.1. Note, the indices for the lower diagonal is (i + 1,
      # i)
      self.target = [10]*self.NN
      self.Qdiag = [1]*self.NN
      self.Qlowerdiag = [0.1]*(self.NN - 1)
 
      super().__init__()
 
   def buildModel(self):
      """
      adding the variables and constraints to the model
      @ingroup PYTHON1THREAD_QP3
      """
      # adding the variables to the model
      for i in range(self.NN):
         self.addVariable(0.0, pyconopt.CONOPT_INF)
 
      # adding the constraints to the model
      # the first constraint is the quadratic objective
      varidx = list(range(self.NN))
      zeros = [0]*self.NN
      ones = [1]*self.NN
      self.addConstraint(pyconopt.ConstraintType_Free, 0.0, varidx, zeros,
                         ones)
 
      # the second constraint is the summation constraint: sum(x) = 1
      self.addConstraint(pyconopt.ConstraintType_Eq, 1.0, varidx, ones,
                    zeros)
 
      # setting the objective constraint
      self.setObjectiveElement(pyconopt.ObjectiveElement_Constraint, 0)
 
      # setting the optimisation direction
      self.setOptimizationSense(pyconopt.Sense_Minimize)
 
      # setting the lagrangian structure for computing the second derivative
      self.setLagrangianStructure()
 
   def setLagrangianStructure(self):
      """
      @ingroup PYTHON1THREAD_QP3
      """
      hessianrow = [f for x in range(self.NN - 1) for f in [x, x + 1]] \
            + [self.NN - 1]
      hessiancol = [f for x in range(self.NN - 1) for f in [x, x]] \
            + [self.NN - 1]
 
      self.setSDLagrangianStructure(hessianrow, hessiancol)
 
   def evaluateNonlinearTerm(self, x, rowno, ignerr, thread):
      """
      @ingroup PYTHON1THREAD_QP3
      """
      # only the objective function is nonlinear, so this function will only be
      # called when rowno == 0
      g = 0
      if rowno == 0:
         g += sum([(x[i] - self.target[i])*q*(x[i] - self.target[i])
                   for i, q in enumerate(self.Qdiag)])
         g += 2*sum([(x[i + 1] - self.target[i + 1])*q*(x[i] - self.target[i])
                   for i, q in enumerate(self.Qlowerdiag)])
 
      return g/2
 
   def evaluateNonlinearJacobian(self, x, rowno, jacnum, ignerr, thread):
      """
      @ingroup PYTHON1THREAD_QP3
      """
      jac = [0]*self.NN
      if rowno == 0:
         for i in range(self.NN):
            jac[i] += self.Qdiag[i]*(x[i] - self.target[i])
            if i < self.NN - 1:
               jac[i + 1] += self.Qlowerdiag[i]*(x[i] - self.target[i])
               jac[i] += self.Qlowerdiag[i]*(x[i + 1] - self.target[i + 1])
 
 
      return jac
 
   def evaluateSDLagrangian(self, x, u, hessianrow, hessiancol):
      """
      @ingroup PYTHON1THREAD_QP3
      """
      hessian = [x*u[0] for i in range(self.NN - 1) \
            for x in [self.Qdiag[i], self.Qlowerdiag[i]]] \
            + [self.Qdiag[-1]*u[0]]
 
      return hessian
 
 
if __name__ == "__main__":
   name = os.path.basename(__file__)[:-3]
 
   conopt = pyconopt.Conopt(name)
   model = QPModelData()
   msghdlr = std.TutMessageHandler(name)
 
   model.buildModel()
 
   conopt.loadModel(model)
   conopt.setMessageHandler(msghdlr)
 
   # getting the license variables
   license_int_1 = os.environ.get('LICENSE_INT_1', None)
   license_int_2 = os.environ.get('LICENSE_INT_2', None)
   license_int_3 = os.environ.get('LICENSE_INT_3', None)
   license_text = os.environ.get('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, 59978.0, coi_error, 0.001)
 
   sys.exit(retcode)