pindyck.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 PDModelData(co.ModelData):
   """ """
 
   def __init__(self):
      self.CONS_SEQ = 0
      self.CONS_DREV = 1
      self.T = 16  # Number of time periods
      self.Vpp = 7  # Variables per period
      self.Epp = 6  # Equations per period
 
      # index arrays for the variables
      self.vartd = []
      self.varcs = []
      self.vars = []
      self.vard = []
      self.varr = []
      self.varp = []
      self.varrev = []
 
      self.consmapping = {}
      super().__init__()
 
   def buildModel(self):
      """
      adding the variables and constraints to the model
      @ingroup PYTHON1THREAD_PINDYCK
      """
 
      demand = [1.0 + 2.3 * pow(1.015, t) for t in range(self.T)]
 
      # adding the variables to the model
      for t in range(self.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)
 
      # adding the objective
      objcoeff = [pow(1.05, 1 - (t + 1)) for t in range(self.T)]
      objnlflag = [0] * self.T
      objidx = self.addConstraint(
         co.ConstraintType_Free, 0.0, self.varrev, objcoeff, objnlflag
      )
 
      # adding the constraints to the model
      for t in range(self.T):
         if t == 0:
            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:
            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],
            )
 
         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.consmapping[considx] = (self.CONS_SEQ, t)
 
         if t == 0:
            self.addConstraint(
               co.ConstraintType_Eq,
               0.0,
               [self.varcs[t], self.vars[t]],
               [1.0, -1.0],
               [0, 0],
            )
         else:
            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.addConstraint(
            co.ConstraintType_Eq,
            0.0,
            [self.vard[t], self.vartd[t], self.vars[t]],
            [1.0, -1.0, 1.0],
            [0, 0, 0],
         )
 
         if t == 0:
            self.addConstraint(
               co.ConstraintType_Eq,
               500.0,
               [self.varr[t], self.vard[t]],
               [1.0, 1.0],
               [0, 0],
            )
         else:
            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],
            )
 
         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.consmapping[considx] = (self.CONS_DREV, t)
 
      # setting the objective constraint
      self.setObjectiveElement(co.ObjectiveElement_Constraint, objidx)
 
      # setting the optimisation direction
      self.setOptimizationSense(co.Sense_Maximize)
 
   def evaluateNonlinearTerm(self, x, rowno, ignerr, thread):
      """
      @copydoc conopt.ModelData.evaluateNonlinearTerm
      @ingroup PYTHON1THREAD_PINDYCK
      """
      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_PINDYCK
      """
      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
 
 
if __name__ == '__main__':
   name = os.path.basename(__file__)[:-3]
 
   conopt = co.Conopt(name)
   model = PDModelData()
   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, 1170.486, coi_error, tol=0.001)
 
   sys.exit(retcode)