Refactor DCOPF for performance

ams/core/matprocessor.py

@@ -51,7 +51,6 @@ def __init__(self,
  v: Optional[np.ndarray] = None,
  sparse: Optional[bool] = False,
  ):
- Param.__init__(self, name=name, info=info)
  self.name = name
  self.tex_name = tex_name if (tex_name is not None) else name
  self.info = info
@@ -107,36 +106,18 @@ def __init__(self, system):
  self.PTDF = MParam(name='PTDF', tex_name=r'P_{TDF}',
  info='Power transfer distribution factor',
  v=None)
- self.pl = MParam(name='pl', tex_name=r'p_l',
- info='Nodal active load',
- v=None, sparse=True)
- self.ql = MParam(name='ql', tex_name=r'q_l',
- info='Nodal reactive load',
- v=None, sparse=True)
  self.Cft = MParam(name='Cft', tex_name=r'C_{ft}',
  info='Connectivity matrix',
  v=None, sparse=True)
- self.Cfti = MParam(name='Cfti', tex_name=r'C_{ft}^i',
- info='Connectivity matrix',
- v=None, sparse=False)
  self.Cg = MParam(name='Cg', tex_name=r'C_g',
  info='Generator connectivity matrix',
  v=None, sparse=True)
- self.Cgi = MParam(name='Cgi', tex_name=r'C_g^i',
- info='Generator connectivity matrix',
- v=None, sparse=True)
  self.Cs = MParam(name='Cs', tex_name=r'C_s',
  info='Slack connectivity matrix',
  v=None, sparse=True)
- self.Csi = MParam(name='Csi', tex_name=r'C_s^i',
- info='Slack connectivity matrix',
- v=None, sparse=True)
  self.Cl = MParam(name='Cl', tex_name=r'Cl',
  info='Load connectivity matrix',
  v=None, sparse=True)
- self.Cli = MParam(name='Cli', tex_name=r'Cl^i',
- info='Load connectivity matrix',
- v=None, sparse=True)
 
  def make(self):
  """
@@ -146,10 +127,8 @@ def make(self):
  """
  system = self.system
  ppc = system2ppc(system)
-
  self.PTDF._v = makePTDF(ppc['baseMVA'], ppc['bus'], ppc['branch'])
  _, _, _, _, self.Cft._v = makeBdc(ppc['baseMVA'], ppc['bus'], ppc['branch'])
- self.Cfti._v = np.linalg.pinv(self.Cft.v)
 
  gen_bus = system.StaticGen.get(src='bus', attr='v',
  idx=system.StaticGen.get_idx())
@@ -158,23 +137,16 @@ def make(self):
  all_bus = system.Bus.idx.v
  load_bus = system.StaticLoad.get(src='bus', attr='v',
  idx=system.StaticLoad.get_idx())
- idx_PD = system.PQ.find_idx(keys="bus", values=all_bus,
- allow_none=True, default=None)
- self.pl._v = c_sparse(system.PQ.get(src='p0', attr='v', idx=idx_PD))
- self.ql._v = c_sparse(system.PQ.get(src='q0', attr='v', idx=idx_PD))
 
  row, col = np.meshgrid(all_bus, slack_bus)
  Cs_v = (row == col).astype(int)
  self.Cs._v = c_sparse(Cs_v)
- self.Csi._v = c_sparse(np.linalg.pinv(Cs_v))
  row, col = np.meshgrid(all_bus, gen_bus)
  Cg_v = (row == col).astype(int)
  self.Cg._v = c_sparse(Cg_v)
- self.Cgi._v = c_sparse(np.linalg.pinv(Cg_v))
  row, col = np.meshgrid(all_bus, load_bus)
  Cl_v = (row == col).astype(int)
  self.Cl._v = c_sparse(Cl_v)
- self.Cli._v = c_sparse(np.linalg.pinv(Cl_v))
 
  return True
 

ams/routines/dcopf.py

@@ -89,29 +89,19 @@ def __init__(self, system, config):
  self.rate_a = RParam(info='long-term flow limit',
  name='rate_a', tex_name=r'R_{ATEA}',
  unit='MVA', model='Line',)
-
  # --- connection matrix ---
- sparse = True
  self.Cg = RParam(info='Gen connection matrix',
  name='Cg', tex_name=r'C_{G}',
  model='mats', src='Cg',
- no_parse=True, sparse=sparse,)
- self.Cgi = RParam(info='Gen connection matrix inverse',
- name='Cgi', tex_name=r'C_{G}^{-1}',
- model='mats', src='Cgi',
- no_parse=True, sparse=sparse,)
+ no_parse=True, sparse=True,)
  self.Cl = RParam(info='Load connection matrix',
  name='Cl', tex_name=r'C_{L}',
  model='mats', src='Cl',
- no_parse=True, sparse=sparse,)
- self.Cli = RParam(info='Load connection matrix inverse',
- name='Cli', tex_name=r'C_{L}^{-1}',
- model='mats', src='Cli',
- no_parse=True, sparse=sparse,)
+ no_parse=True, sparse=True,)
  self.Cft = RParam(info='Line connection matrix',
  name='Cft', tex_name=r'C_{ft}',
  model='mats', src='Cft',
- no_parse=True, sparse=sparse,)
+ no_parse=True, sparse=True,)
  self.PTDF = RParam(info='Power Transfer Distribution Factor',
  name='PTDF', tex_name=r'P_{TDF}',
  model='mats', src='PTDF',
@@ -123,23 +113,29 @@ def __init__(self, system, config):
  name='pg', tex_name=r'p_{G}',
  model='StaticGen', src='p',
  v0=self.pg0)
+ # NOTE: Var bounds need to set separately
  pglb = '-pg + mul(nctrle, pg0) + mul(ctrle, pmin)'
  self.pglb = Constraint(name='pglb', info='pg min',
  e_str=pglb, type='uq',)
  pgub = 'pg - mul(nctrle, pg0) - mul(ctrle, pmax)'
  self.pgub = Constraint(name='pgub', info='pg max',
  e_str=pgub, type='uq',)
  # --- bus ---
- self.aBus = Var(info='Bus voltage angle',
- name='aBus', tex_name=r'\theta_{Bus}',
- unit='rad',
- model='Bus', src='a',)
- self.vBus = Var(info='Bus voltage magnitude',
- name='vBus', tex_name=r'v_{Bus}',
- unit='p.u.',
- model='Bus', src='v')
+ self.png = Var(info='Bus active power from gen',
+ unit='p.u.',
+ name='pg', tex_name=r'p_{Bus,G}',
+ model='Bus',)
+ self.pnd = Var(info='Bus active power from load',
+ unit='p.u.',
+ name='pd', tex_name=r'p_{Bus,D}',
+ model='Bus',)
+ self.pngb = Constraint(name='pngb', type='eq',
+ e_str='Cg@png - pg',
+ info='Bus active power from gen',)
+ self.pndb = Constraint(name='pndb', type='eq',
+ e_str='Cl@pnd - pd',
+ info='Bus active power from load',)
  # --- line ---
- # NOTE: Var bounds need to set separately
  # NOTE: `ug*pmin` results in unexpected error
  self.plf = Var(info='Line active power',
  name='plf', tex_name=r'p_{L}',
@@ -154,7 +150,7 @@ def __init__(self, system, config):
  type='eq',)
  self.pnb = Constraint(name='pnb', type='eq',
  info='nodal power injection',
- e_str='PTDF@(Cgi@pg - Cli@pd) - plf',)
+ e_str='PTDF@(png - pnd) - plf',)
  # --- objective ---
  obj = 'sum(mul(c2, power(pg, 2)))' \
  '+ sum(mul(c1, pg))' \
@@ -167,17 +163,7 @@ def solve(self, **kwargs):
  """
  Solve the routine optimization model.
  """
- res = self.om.prob.solve(**kwargs)
- # str output maens solver failed
- if isinstance(res, str):
- return res
- # estimate aBus
- aBus = np.matmul(self.system.mats.Cfti.v,
- np.multiply(self.x.v, self.plf.v))
- self.aBus.optz.value = aBus
- # set vBus to 1
- self.vBus.optz.value = np.ones(self.vBus.shape)
- return res
+ return self.om.prob.solve(**kwargs)
 
  def run(self, no_code=True, **kwargs):
  """

ams/routines/dopf.py

@@ -26,10 +26,19 @@ def __init__(self, system, config):
  self.info = 'Linearzied distribution OPF'
  self.type = 'DED'
 
- # --- params ---
+ # -- Data Section --
+ # --- generator ---
+ self.qmax = RParam(info='generator maximum reactive power',
+ name='qmax', tex_name=r'q_{max}', unit='p.u.',
+ model='StaticGen', src='qmax',)
+ self.qmin = RParam(info='generator minimum reactive power',
+ name='qmin', tex_name=r'q_{min}', unit='p.u.',
+ model='StaticGen', src='qmin',)
+ # --- load ---
  self.qd = RParam(info='reactive demand',
  name='qd', tex_name=r'q_{d}', unit='p.u.',
  model='StaticLoad', src='q0',)
+ # --- bus ---
  self.vmax = RParam(info="Bus voltage upper limit",
  name='vmax', tex_name=r'\overline{v}',
  unit='p.u.',
@@ -39,16 +48,12 @@ def __init__(self, system, config):
  name='vmin', tex_name=r'\underline{v}',
  unit='p.u.',
  model='Bus', src='vmin', no_parse=True,)
+ # --- line ---
  self.r = RParam(info='line resistance',
  name='r', tex_name='r', unit='p.u.',
  model='Line', src='r')
- self.qmax = RParam(info='generator maximum reactive power',
- name='qmax', tex_name=r'q_{max}', unit='p.u.',
- model='StaticGen', src='qmax',)
- self.qmin = RParam(info='generator minimum reactive power',
- name='qmin', tex_name=r'q_{min}', unit='p.u.',
- model='StaticGen', src='qmin',)
- # --- vars ---
+ # --- Model Section ---
+ # --- generator ---
  self.qg = Var(info='Gen reactive power',
  name='qg', tex_name=r'q_{g}', unit='p.u.',
  model='StaticGen', src='q',)
@@ -58,7 +63,7 @@ def __init__(self, system, config):
  self.qgub = Constraint(name='qgub', type='uq',
  info='qg max',
  e_str='qg - mul(ug, qmax)',)
-
+ # --- bus ---
  self.vsq = Var(info='square of Bus voltage',
  name='vsq', tex_name=r'v^{2}', unit='p.u.',
  model='Bus',)
@@ -70,21 +75,17 @@ def __init__(self, system, config):
  info='Voltage lower limit',
  e_str='-vsq + vmin**2',
  type='uq',)
-
+ # --- line ---
  self.qlf = Var(info='line reactive power',
  name='qlf', tex_name=r'q_{lf}',
  unit='p.u.', model='Line',)
-
- # --- constraints ---
- self.pnb.e_str = 'PTDF@(Cgi@pg - Cli@pd) - plf'
-
- self.qb = Constraint(name='qb', info='reactive power balance',
- e_str='sum(qd) - sum(qg)',
- type='eq',)
- self.lvd = Constraint(name='lvd',
- info='line voltage drop',
- e_str='Cft@vsq - (r * plf + x * qlf)',
- type='eq',)
+ self.lvd = Constraint(info='line voltage drop',
+ name='lvd', type='eq',
+ e_str='Cft@vsq - (r * plf + x * qlf)',)
+ # --- power balance ---
+ self.qb = Constraint(info='reactive power balance',
+ name='qb', type='eq',
+ e_str='sum(qd) - sum(qg)',)
 
  # --- objective ---
  # NOTE: no need to revise objective function

ams/routines/ed.py

@@ -64,6 +64,7 @@ def __init__(self) -> None:
  name='sd', tex_name=r's_{d}',
  src='sd', model='EDTSlot')
 
+ # NOTE: update timeslot.model in dispatch model if necessary
  self.timeslot = RParam(info='Time slot for multi-period ED',
  name='timeslot', tex_name=r't_{s,idx}',
  src='idx', model='EDTSlot',
@@ -104,6 +105,10 @@ def __init__(self) -> None:
  self.rate_a.expand_dims = 1
  self.x.expand_dims = 1
 
+ # NOTE: extend pg to 2D matrix: row for gen and col for timeslot
+ self.pg.horizon = self.timeslot
+ self.pg.info = '2D Gen power'
+
 class ED(RTED):
  """
  DC-based multi-period economic dispatch (ED).
@@ -147,16 +152,14 @@ def __init__(self, system, config):
  self.dud.expand_dims = 1
  self.ddd.expand_dims = 1
 
- # --- params ---
+ # --- Data Section ---
  self.ugt = NumOp(u=self.ug, fun=np.transpose,
  name='ugt', tex_name=r'u_{g}',
  info='input ug transpose',
  no_parse=True)
 
- # --- vars ---
- # NOTE: extend pg to 2D matrix, where row is gen and col is timeslot
- self.pg.horizon = self.timeslot
- self.pg.info = '2D Gen power'
+ # --- Model Section ---
+ # --- gen ---
  self.ctrle.u2 = self.ugt
  self.nctrle.u2 = self.ugt
  pglb = '-pg + mul(mul(nctrle, pg0), tlv) '
@@ -166,34 +169,31 @@ def __init__(self, system, config):
  pgub += '- mul(mul(ctrle, tlv), pmax)'
  self.pgub.e_str = pgub
 
- self.plf.horizon = self.timeslot
- self.plf.info = '2D Line flow'
-
  self.pru.horizon = self.timeslot
  self.pru.info = '2D RegUp power'
-
  self.prd.horizon = self.timeslot
  self.prd.info = '2D RegDn power'
 
  self.prs.horizon = self.timeslot
-
- self.aBus.horizon = self.timeslot
- self.aBus.info = '2D Bus angle'
-
- # --- constraints ---
- # NOTE: Spg @ pg returns a row vector
- self.pb.e_str = '- gs @ pg + pdsz' # power balance
-
  self.prsb.e_str = 'mul(ugt, mul(pmax, tlv) - pg) - prs'
  self.rsr.e_str = '-gs@prs + dsr'
 
- # --- bus power injection ---
- self.pnb.e_str = 'PTDF@(Cgi@pg - Cli@pds) - plf'
-
- # --- line limits ---
+ # --- line ---
+ self.plf.horizon = self.timeslot
+ self.plf.info = '2D Line flow'
  self.plflb.e_str = '-plf - mul(rate_a, tlv)'
  self.plfub.e_str = 'plf - mul(rate_a, tlv)'
 
+ # --- power balance ---
+ # NOTE: Spg @ pg returns a row vector
+ self.pb.e_str = '- gs @ pg + pdsz'
+
+ self.png.horizon = self.timeslot
+ self.pnb.e_str = 'PTDF@(png - pnd) - plf'
+
+ self.pnd.horizon = self.timeslot
+ self.pndb.e_str = 'Cl@pnd - pds'
+
  # --- ramping ---
  self.rbu.e_str = 'gs@mul(ugt, pru) - mul(dud, tlv)'
  self.rbd.e_str = 'gs@mul(ugt, prd) - mul(ddd, tlv)'
@@ -228,23 +228,6 @@ def dc2ac(self, **kwargs):
  """
  return NotImplementedError
 
- def solve(self, **kwargs):
- """
- Solve the routine optimization model.
- """
- res = self.om.prob.solve(**kwargs)
- # str output maens solver failed
- if isinstance(res, str):
- return res
- # estimate aBus
- x_tlv = np.matmul(self.x.v, self.tlv.v)
- aBus = np.matmul(self.system.mats.Cfti.v,
- np.multiply(x_tlv, self.plf.v))
- self.aBus.optz.value = aBus
- # set vBus to 1
- self.vBus.optz.value = np.ones(self.vBus.shape)
- return None
-
  def unpack(self, **kwargs):
  """
  Multi-period dispatch will not unpack results from

ams/routines/routine.py

@@ -298,8 +298,10 @@ def init(self, force=True, no_code=True, **kwargs):
  msg = f"{self.class_name} data check failed, setup may run into error!"
  logger.warning(msg)
  self._constr_check()
- # FIXME: build the system matrices every init might slow down
+ t_mat, _ = elapsed()
  self.system.mats.make()
+ _, s_mat = elapsed(t_mat)
+ logger.debug(f"Built system matrices in {s_mat}.")
  t_setup, _ = elapsed()
  results = self.om.setup(no_code=no_code)
  _, s_setup = elapsed(t_setup)