qp3.py

Go to the documentation of this file.

 
import os
import sys
 
import conopt as co
 
sys.path.append('../common/')
import std
 
 
class QPModelData(co.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, co.Conopt.Infinity)
 
      # 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(co.ConstraintType_Free, 0.0, varidx, zeros, ones)
 
      # the second constraint is the summation constraint: sum(x) = 1
      self.addConstraint(co.ConstraintType_Eq, 1.0, varidx, ones, zeros)
 
      # setting the objective constraint
      self.setObjectiveElement(co.ObjectiveElement_Constraint, 0)
 
      # setting the optimisation direction
      self.setOptimizationSense(co.Sense_Minimize)
 
      # setting the lagrangian structure for computing the second derivative
      self.setLagrangianStructure()
 
   def setLagrangianStructure(self):
      """
      define the lagrangian structure and supplies it to CONOPT
 
      @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):
      """
      @copydoc conopt.ModelData.evaluateNonlinearTerm
      @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):
      """
      @copydoc conopt.ModelData.evaluateNonlinearJacobian
      @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):
      """
      @copydoc conopt.ModelData.evaluateSDLagrangian
      @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 = co.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('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, 59978.0, coi_error, 0.001)
 
   sys.exit(retcode)