pinadd2ddir.py

Go to the documentation of this file.

 
import os
import sys
import math
 
import conopt as co
 
sys.path.append('../common/')
import std
 
 
class PinAdd2DDirModelData(co.ModelData):
   """ """
 
   def __init__(self):
      self.CONS_SEQ = 0
      self.CONS_DREV = 1
 
      self.T = 0
 
      # index arrays for the variables
      self.vartd = []
      self.varcs = []
      self.vars = []
      self.vard = []
      self.varr = []
      self.varp = []
      self.varrev = []
 
      # index arrays for the constraints
      self.consttdeq = []
      self.constseq = []
      self.constcseq = []
      self.constdeq = []
      self.constreq = []
      self.constdrev = []
 
      self.consmapping = {}
 
      self.hessianstart = []
      super().__init__()
 
   def buildModel(self, T, xkeep, xstat, estat):
      """
      adding the variables and constraints to the model
      @ingroup PYTHON1THREAD_PINADD2DDIR
      """
 
      self.T = T
      demand = [1.0 + 2.3 * pow(1.015, t) for t in range(T)]
 
      # adding the variables to the model
      for t in range(T):
         # td:  Lower=0, Upper=inf, Curr=18
         varidx = self.addVariable(0, co.Conopt.Infinity, 18)
         self.vartd.append(varidx)
 
         # cs:  Lower=0, Upper=inf, Curr=7*t
         varidx = self.addVariable(0, co.Conopt.Infinity, 7 * (t + 1))
         self.varcs.append(varidx)
 
         # s:   Lower=0, Upper=inf, Curr=7
         varidx = self.addVariable(0, co.Conopt.Infinity, 7)
         self.vars.append(varidx)
 
         # d:   Lower=0, Upper=inf, Curr=td-s
         varidx = self.addVariable(
            0,
            co.Conopt.Infinity,
            self.getVariable(self.vartd[t]).curr
            - self.getVariable(self.vars[t]).curr,
         )
         self.vard.append(varidx)
 
         # r:   Lower=1, Upper=inf, Curr=r(t-1)-d
         if t > 0:
            varidx = self.addVariable(
               1,
               co.Conopt.Infinity,
               self.getVariable(self.varr[t - 1]).curr
               - self.getVariable(self.vard[t]).curr,
            )
         else:
            varidx = self.addVariable(
               1,
               co.Conopt.Infinity,
               500 - self.getVariable(self.vard[t]).curr,
            )
         self.varr.append(varidx)
 
         # p:   Lower=1, Upper=inf, Curr=14
         varidx = self.addVariable(1, co.Conopt.Infinity, 14)
         self.varp.append(varidx)
 
         # rev: Lower=-inf, Upper=inf, Curr=0
         varidx = self.addVariable(-co.Conopt.Infinity, co.Conopt.Infinity)
         self.varrev.append(varidx)
 
      # This is a restart: Use the initial values from last solve for
      # the variables in the first periods and extrapolate the last
      # period using a linear extrapolation between the last two
      # periods.
      if len(xkeep) > 0:
         for i in range(len(xkeep)):
            self.getVariable(i).curr = xkeep[i]
 
         # linear extrapolation for the last period of variable td
         self.getVariable(self.vartd[T - 1]).curr = (
            2 * self.getVariable(self.vartd[T - 2]).curr
            - self.getVariable(self.vartd[T - 3]).curr
         )
         # cs
         self.getVariable(self.varcs[T - 1]).curr = (
            2 * self.getVariable(self.varcs[T - 2]).curr
            - self.getVariable(self.varcs[T - 3]).curr
         )
         # s
         self.getVariable(self.vars[T - 1]).curr = (
            2 * self.getVariable(self.vars[T - 2]).curr
            - self.getVariable(self.vars[T - 3]).curr
         )
         # d
         self.getVariable(self.vard[T - 1]).curr = (
            2 * self.getVariable(self.vard[T - 2]).curr
            - self.getVariable(self.vard[T - 3]).curr
         )
         # r
         self.getVariable(self.varr[T - 1]).curr = (
            2 * self.getVariable(self.varr[T - 2]).curr
            - self.getVariable(self.varr[T - 3]).curr
         )
         # p
         self.getVariable(self.varp[T - 1]).curr = (
            2 * self.getVariable(self.varp[T - 2]).curr
            - self.getVariable(self.varp[T - 3]).curr
         )
         # rev
         self.getVariable(self.varrev[T - 1]).curr = (
            2 * self.getVariable(self.varrev[T - 2]).curr
            - self.getVariable(self.varrev[T - 3]).curr
         )
 
         # The variables from the last solve are given the old status and
         # the variables in the new period are given those in the last
         # period. Similarly with the Equation status:
         for i in range(len(xstat)):
            self.getVariable(i).varstatus = xstat[i]
 
         # td
         self.getVariable(self.vartd[T - 1]).varstatus = self.getVariable(
            self.vartd[T - 2]
         ).varstatus
         # cs
         self.getVariable(self.varcs[T - 1]).varstatus = self.getVariable(
            self.varcs[T - 2]
         ).varstatus
         # s
         self.getVariable(self.vars[T - 1]).varstatus = self.getVariable(
            self.vars[T - 2]
         ).varstatus
         # d
         self.getVariable(self.vard[T - 1]).varstatus = self.getVariable(
            self.vard[T - 2]
         ).varstatus
         # r
         self.getVariable(self.varr[T - 1]).varstatus = self.getVariable(
            self.varr[T - 2]
         ).varstatus
         # p
         self.getVariable(self.varp[T - 1]).varstatus = self.getVariable(
            self.varp[T - 2]
         ).varstatus
         # rev
         self.getVariable(self.varrev[T - 1]).varstatus = self.getVariable(
            self.varrev[T - 2]
         ).varstatus
 
      # adding the objective
      objcoeff = [pow(1.05, 1 - (t + 1)) for t in range(T)]
      objnlflag = [0] * T
      objidx = self.addConstraint(
         co.ConstraintType_Free, 0.0, self.varrev, objcoeff, objnlflag
      )
 
      # adding the constraints to the model
      for t in range(T):
         if t == 0:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               demand[t] + 0.87 * 18.0,
               [self.vartd[t], self.varp[t]],  # the variables
               [1.0, 0.13],  # the coefficients
               [0, 0],  # nlflags
            )
         else:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               demand[t],
               [self.vartd[t], self.vartd[t - 1], self.varp[t]],
               [1.0, -0.87, 0.13],
               [0, 0, 0],
            )
         self.consttdeq.append(considx)
 
         if t == 0:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               0.75 * 6.5,
               [self.vars[t], self.varp[t], self.varcs[t]],
               [1.0, 0.0, 0.0],
               [0, 1, 1],
            )
         else:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               0.0,
               [
                  self.vars[t],
                  self.vars[t - 1],
                  self.varp[t],
                  self.varcs[t],
               ],
               [1.0, -0.75, 0.0, 0.0],
               [0, 0, 1, 1],
            )
         self.constseq.append(considx)
         self.consmapping[considx] = (self.CONS_SEQ, t)
 
         if t == 0:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               0.0,
               [self.varcs[t], self.vars[t]],
               [1.0, -1.0],
               [0, 0],
            )
         else:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               0.0,
               [self.varcs[t], self.varcs[t - 1], self.vars[t]],
               [1.0, -1.0, -1.0],
               [0, 0, 0],
            )
         self.constcseq.append(considx)
 
         considx = self.addConstraint(
            co.ConstraintType_Eq,
            0.0,
            [self.vard[t], self.vartd[t], self.vars[t]],
            [1.0, -1.0, 1.0],
            [0, 0, 0],
         )
         self.constdeq.append(considx)
 
         if t == 0:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               500.0,
               [self.varr[t], self.vard[t]],
               [1.0, 1.0],
               [0, 0],
            )
         else:
            considx = self.addConstraint(
               co.ConstraintType_Eq,
               0.0,
               [self.varr[t], self.varr[t - 1], self.vard[t]],
               [1.0, -1.0, 1.0],
               [0, 0, 0],
            )
         self.constreq.append(considx)
 
         considx = self.addConstraint(
            co.ConstraintType_Eq,
            0.0,
            [self.varrev[t], self.vard[t], self.varp[t], self.varr[t]],
            [1.0, 0.0, 0.0, 0.0],
            [0, 1, 1, 1],
         )
         self.constdrev.append(considx)
         self.consmapping[considx] = (self.CONS_DREV, t)
 
      # Restore values from previous runs
      if len(estat) > 0:
         # the equation status
         for i in range(len(estat)):
            self.getConstraint(i).slackstatus = estat[i]
 
         # tdeq
         self.getConstraint(
            self.consttdeq[T - 1]
         ).slackstatus = self.getConstraint(self.consttdeq[T - 2]).slackstatus
         # seq
         self.getConstraint(
            self.constseq[T - 1]
         ).slackstatus = self.getConstraint(self.constseq[T - 2]).slackstatus
         # cseq
         self.getConstraint(
            self.constcseq[T - 1]
         ).slackstatus = self.getConstraint(self.constcseq[T - 2]).slackstatus
         # deq
         self.getConstraint(
            self.constdeq[T - 1]
         ).slackstatus = self.getConstraint(self.constdeq[T - 2]).slackstatus
         # req
         self.getConstraint(
            self.constreq[T - 1]
         ).slackstatus = self.getConstraint(self.constreq[T - 2]).slackstatus
         # drev
         self.getConstraint(
            self.constdrev[T - 1]
         ).slackstatus = self.getConstraint(self.constdrev[T - 2]).slackstatus
 
      # setting the objective constraint
      self.setObjectiveElement(co.ObjectiveElement_Constraint, objidx)
 
      # setting the optimisation direction
      self.setOptimizationSense(co.Sense_Maximize)
 
      hsrow = []
      hscol = []
 
      for t in range(T):
         # storing the start index for the hessian in each time period.
         # This will provide a mapping that can be used in the hessian
         # evaluation
         self.hessianstart.append(len(hscol))
 
         hscol.append(self.varcs[t])
         hsrow.append(self.varcs[t])
         hscol.append(self.varcs[t])
         hsrow.append(self.varp[t])
         hscol.append(self.vard[t])
         hsrow.append(self.varr[t])
         hscol.append(self.vard[t])
         hsrow.append(self.varp[t])
         hscol.append(self.varr[t])
         hsrow.append(self.varr[t])
      self.setSDLagrangianStructure(hsrow, hscol)
 
   def evaluateNonlinearTerm(self, x, rowno, ignerr, thread):
      """
      @copydoc conopt.ModelData.evaluateNonlinearTerm
      @ingroup PYTHON1THREAD_PINADD2DDIR
      """
      if rowno not in self.consmapping:
         return
 
      # unpack consset and time index
      consset, t = self.consmapping[rowno]
      assert consset in (self.CONS_SEQ, self.CONS_DREV)
 
      g = 0
      if consset == self.CONS_SEQ:
         h1 = 1.1 + 0.1 * x[self.varp[t]]
         h2 = pow(1.02, -x[self.varcs[t]] / 7.0)
         g = -h1 * h2
 
      elif consset == self.CONS_DREV:
         g = -x[self.vard[t]] * (x[self.varp[t]] - 250.0 / x[self.varr[t]])
 
      return g
 
   def evaluateNonlinearJacobian(self, x, rowno, jacnum, ignerr, thread):
      """
      @copydoc conopt.ModelData.evaluateNonlinearJacobian
      @ingroup PYTHON1THREAD_PINADD2DDIR
      """
      if rowno not in self.consmapping:
         return
 
      # unpack consset and time index
      consset, t = self.consmapping[rowno]
      assert consset in (self.CONS_SEQ, self.CONS_DREV)
 
      jac = []
      if consset == self.CONS_SEQ:
         h1 = 1.1 + 0.1 * x[self.varp[t]]
         h2 = pow(1.02, -x[self.varcs[t]] / 7.0)
         jac.append(h1 * h2 * math.log(1.02) / 7.0)
         jac.append(-h2 * 0.1)
 
      elif consset == self.CONS_DREV:
         jac.append(-(x[self.varp[t]] - 250.0 / x[self.varr[t]]))
         jac.append(-x[self.vard[t]] * 250.0 / pow(x[self.varr[t]], 2))
         jac.append(-x[self.vard[t]])
 
      return jac
 
   def evaluateDirectionalSD(self, x, dx, rowno, jacnum, thread):
      """
      @copydoc conopt.ModelData.evaluateDirectionalSD
      @ingroup PYTHON1THREAD_PINADD2DDIR
      """
      # returning an iterator based on the row number.
      # NOTE: a map is just an array of pairs. So for the iterator, the
      # first element is the key and the second is the value. In our case,
      # then value is also a pair.
      if rowno not in self.consmapping:
         return
 
      # storing the value (a pair) as a separate variable to make the
      # extraction easier
      # unpack consset and time index
      consset, t = self.consmapping[rowno]
      assert consset in (self.CONS_SEQ, self.CONS_DREV)
      dirsd = [0 for i in range(jacnum)]
      if consset == self.CONS_SEQ:
         # seq equation. Nonlinear term = -(1.1+0.1*p)*1.02**(-cs/7)
         # Hessian (diagonal and lower triangle): A total of 2 nonzero elements per period
         #                         1                                  5
         #                         cs                                  p
         # 1 : cs  -(1.1+0.1p)*1.02**(-cs/7)*(log(1.02)/7)**2
         # 5 : p           0.1*1.02**(-cs/7)*(log(1.02)72)
 
         cs = x[self.varcs[t]]
         p = x[self.varp[t]]
         h1 = 1.1 + 0.1 * p
         h2 = pow(1.02, -cs / 7.0)
         h11 = -h1 * h2 * pow(math.log(1.02) / 7.0, 2)
         h15 = 0.1 * h2 * (math.log(1.02) / 7.0)
         dirsd[self.varcs[t]] = h11 * dx[self.varcs[t]] + h15 * dx[self.varp[t]]
         dirsd[self.varp[t]] = h15 * dx[self.varcs[t]]
 
      elif consset == self.CONS_DREV:
         # drev equation. Nonlinear term = -d*(p-250/r)
         # Hessian (diagonal and lower triangle): A total of 3 nonzero elements per period
         #                 3           4           5
         #                 d           r           p
         # 3 : d
         # 4 : r     -250/r**2   500*d/r**3
         # 5 : p        -1
         d = x[self.vard[t]]
         r = x[self.varr[t]]
         h34 = -250.0 / pow(r, 2)
         h35 = -1.0
         h44 = 500.0 * d / pow(r, 3)
         dirsd[self.vard[t]] = h34 * dx[self.varr[t]] + h35 * dx[self.varp[t]]
         dirsd[self.varr[t]] = h34 * dx[self.vard[t]] + h44 * dx[self.varr[t]]
         dirsd[self.varp[t]] = h35 * dx[self.vard[t]]
 
      else:
         # Error - this equation is not nonlinear
         print('\nError. Pin_2DDir called with rowno = %d.\n\n', rowno)
         return
 
      return dirsd
 
   def evaluateSDLagrangian(self, x, u, hessianrow, hessiancol):
      """
      @copydoc conopt.ModelData.evaluateSDLagrangian
      @ingroup PYTHON1THREAD_PINADD2DDIR
      """
      numhessian = len(hessianrow)
      hessian = [0 for i in range(numhessian)]
 
      # Normal Evaluation mode
      for t in range(self.T):
         hessidx = self.hessianstart[t]
         cs = x[self.varcs[t]]
         d = x[self.vard[t]]
         r = x[self.varr[t]]
         p = x[self.varp[t]]
 
         h1 = 1.1 + 0.1 * p
         h2 = pow(1.02, -cs / 7.0)
 
         # the second derivative of the SEQ equations
         hessian[hessidx] = (
            -h1 * h2 * pow(math.log(1.02) / 7.0, 2) * u[self.constseq[t]]
         )
         hessian[hessidx + 1] = (
            0.1 * h2 * (math.log(1.02) / 7.0) * u[self.constseq[t]]
         )
 
         # the second derivative of the DREV equations
         hessian[hessidx + 2] = -250.0 / pow(r, 2) * u[self.constdrev[t]]
         hessian[hessidx + 3] = -1.0 * u[self.constdrev[t]]
         hessian[hessidx + 4] = 500.0 * d / pow(r, 3) * u[self.constdrev[t]]
 
      return hessian
 
 
if __name__ == '__main__':
   name = os.path.basename(__file__)[:-3]
 
   conopt = co.Conopt(name)
   msghdlr = std.TutMessageHandler(name)
   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,
      )
 
   expected_sols = [1170.49, 1225.86, 1280.28, 1333.74, 1386.24]
   for i in range(16, 21):
      model = PinAdd2DDirModelData()
 
      xkeep = conopt.getVariableValues()
      xstat = conopt.getVariableStatus()
      estat = conopt.getConstraintStatus()
 
      # building the model
      model.buildModel(
         i, xkeep, xstat, estat
      )  # Add T, and other vectors as parameters
 
      conopt.loadModel(model)
 
      coi_error = conopt.solve()
 
      retcode = std.checkSolve(
         conopt, expected_sols[i - 16], coi_error, tol=0.01
      )
 
   sys.exit(retcode)