pylon/opf.py

   1 #------------------------------------------------------------------------------
   2 # Copyright (C) 1996-2010 Power System Engineering Research Center (PSERC)
   3 # Copyright (C) 2007-2010 Richard Lincoln
   4 #
   5 # Licensed under the Apache License, Version 2.0 (the "License");
   6 # you may not use this file except in compliance with the License.
   7 # You may obtain a copy of the License at
   8 #
   9 #     http://www.apache.org/licenses/LICENSE-2.0
  10 #
  11 # Unless required by applicable law or agreed to in writing, software
  12 # distributed under the License is distributed on an "AS IS" BASIS,
  13 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14 # See the License for the specific language governing permissions and
  15 # limitations under the License.
  16 #------------------------------------------------------------------------------
  17
  18 """ Defines a generalised OPF solver and an OPF model.
  19
  20 Based on opf.m from MATPOWER by Ray Zimmerman, developed at PSERC Cornell.
  21 See U{http://www.pserc.cornell.edu/matpower/} for more info.
  22 """
  23
  24 #------------------------------------------------------------------------------
  25 #  Imports:
  26 #------------------------------------------------------------------------------
  27
  28 import logging
  29 import random
  30
  31 from time import time
  32
  33 from numpy import \
  34     array, pi, diff, Inf, ones, r_, float64, zeros, arctan2, sin, cos
  35
  36 from scipy.sparse import lil_matrix, csr_matrix, hstack
  37
  38 from util import _Named, fair_max
  39 from case import REFERENCE
  40 from generator import PW_LINEAR
  41 from solver import DCOPFSolver, PIPSSolver
  42
  43 #------------------------------------------------------------------------------
  44 #  Logging:
  45 #------------------------------------------------------------------------------
  46
  47 logger = logging.getLogger(__name__)
  48
  49 #------------------------------------------------------------------------------
  50 #  "OPF" class:
  51 #------------------------------------------------------------------------------
  52
  53 class OPF(object):
  54     """ Defines a generalised OPF solver.
  55
  56     Based on opf.m from MATPOWER by Ray Zimmerman, developed at PSERC
  57     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
  58     """
  59
  60     def __init__(self, case, dc=True, ignore_ang_lim=True, opt=None):
  61         """ Initialises a new OPF instance.
  62         """
  63         #: Case under optimisation.
  64         self.case = case
  65
  66         #: Use DC power flow formulation.
  67         self.dc = dc
  68
  69         #: Ignore angle difference limits for branches even if specified.
  70         self.ignore_ang_lim = ignore_ang_lim
  71
  72         #: Solver options (See pips.py for futher details).
  73         self.opt = {} if opt is None else opt
  74
  75     #--------------------------------------------------------------------------
  76     #  Public interface:
  77     #--------------------------------------------------------------------------
  78
  79     def solve(self, solver_klass=None):
  80         """ Solves an optimal power flow and returns a results dictionary.
  81         """
  82         # Start the clock.
  83         t0 = time()
  84
  85         # Build an OPF model with variables and constraints.
  86         om = self._construct_opf_model(self.case)
  87         if om is None:
  88             return {"converged": False, "output": {"message": "No Ref Bus."}}
  89
  90         # Call the specific solver.
  91 #        if self.opt["verbose"]:
  92 #            print '\nPYLON Version %s, %s', "0.4.2", "April 2010"
  93         if solver_klass is not None:
  94             result = solver_klass(om, opt=self.opt).solve()
  95         elif self.dc:
  96 #            if self.opt["verbose"]:
  97 #                print ' -- DC Optimal Power Flow\n'
  98             result = DCOPFSolver(om, opt=self.opt).solve()
  99         else:
 100 #            if self.opt["verbose"]:
 101 #                print ' -- AC Optimal Power Flow\n'
 102             result = PIPSSolver(om, opt=self.opt).solve()
 103
 104         result["elapsed"] = time() - t0
 105
 106         if self.opt.has_key("verbose"):
 107             if self.opt["verbose"]:
 108                 logger.info("OPF completed in %.3fs." % result["elapsed"])
 109
 110         return result
 111
 112     #--------------------------------------------------------------------------
 113     #  Private interface:
 114     #--------------------------------------------------------------------------
 115
 116     def _construct_opf_model(self, case):
 117         """ Returns an OPF model.
 118         """
 119         # Zero the case result attributes.
 120         self.case.reset()
 121
 122         base_mva = case.base_mva
 123
 124         # Check for one reference bus.
 125         oneref, refs = self._ref_check(case)
 126         if not oneref: #return {"status": "error"}
 127             None
 128
 129         # Remove isolated components.
 130         bs, ln, gn = self._remove_isolated(case)
 131
 132         # Update bus indexes.
 133         self.case.index_buses(bs)
 134
 135         # Convert single-block piecewise-linear costs into linear polynomial.
 136         gn = self._pwl1_to_poly(gn)
 137
 138         # Set-up initial problem variables.
 139         Va = self._get_voltage_angle_var(refs, bs)
 140         Pg = self._get_pgen_var(gn, base_mva)
 141
 142         if self.dc: # DC model.
 143             # Get the susceptance matrices and phase shift injection vectors.
 144             B, Bf, Pbusinj, Pfinj = self.case.makeBdc(bs, ln)
 145
 146             # Power mismatch constraints (B*Va + Pg = Pd).
 147             Pmis = self._power_mismatch_dc(bs, gn, B, Pbusinj, base_mva)
 148
 149             # Branch flow limit constraints.
 150             Pf, Pt = self._branch_flow_dc(ln, Bf, Pfinj, base_mva)
 151         else:
 152             # Set-up additional AC-OPF problem variables.
 153             Vm = self._get_voltage_magnitude_var(bs, gn)
 154             Qg = self._get_qgen_var(gn, base_mva)
 155
 156             Pmis, Qmis, Sf, St = self._nln_constraints(len(bs), len(ln))
 157
 158             vl = self._const_pf_constraints(gn, base_mva)
 159
 160             # TODO: Generator PQ capability curve constraints.
 161 #            PQh, PQl = self._pq_capability_curve_constraints(gn)
 162
 163         # Branch voltage angle difference limits.
 164         ang = self._voltage_angle_diff_limit(bs, ln)
 165
 166         if self.dc:
 167             vars = [Va, Pg]
 168             constraints = [Pmis, Pf, Pt, ang]
 169         else:
 170             vars = [Va, Vm, Pg, Qg]
 171             constraints = [Pmis, Qmis, Sf, St, #PQh, PQL,
 172                            vl, ang]
 173
 174         # Piece-wise linear generator cost constraints.
 175         y, ycon = self._pwl_gen_costs(gn, base_mva)
 176
 177         if ycon is not None:
 178             vars.append(y)
 179             constraints.append(ycon)
 180
 181         # Add variables and constraints to the OPF model object.
 182         opf = OPFModel(case)
 183         opf.add_vars(vars)
 184         opf.add_constraints(constraints)
 185
 186         if self.dc: # user data
 187             opf._Bf = Bf
 188             opf._Pfinj = Pfinj
 189
 190         return opf
 191
 192
 193     def _ref_check(self, case):
 194         """ Checks that there is only one reference bus.
 195         """
 196         refs = [bus._i for bus in case.buses if bus.type == REFERENCE]
 197
 198         if len(refs) == 1:
 199             return True, refs
 200         else:
 201             logger.error("OPF requires a single reference bus.")
 202             return False, refs
 203
 204
 205     def _remove_isolated(self, case):
 206         """ Returns non-isolated case components.
 207         """
 208 #        case.deactivate_isolated()
 209         buses = case.connected_buses
 210         branches = case.online_branches
 211         gens = case.online_generators
 212
 213         return buses, branches, gens
 214
 215
 216     def _pwl1_to_poly(self, generators):
 217         """ Converts single-block piecewise-linear costs into linear
 218         polynomial.
 219         """
 220         for g in generators:
 221             if (g.pcost_model == PW_LINEAR) and (len(g.p_cost) == 2):
 222                 g.pwl_to_poly()
 223
 224         return generators
 225
 226     #--------------------------------------------------------------------------
 227     #  Optimisation variables:
 228     #--------------------------------------------------------------------------
 229
 230     def _get_voltage_angle_var(self, refs, buses):
 231         """ Returns the voltage angle variable set.
 232         """
 233         Va = array([b.v_angle * (pi / 180.0) for b in buses])
 234
 235         Vau = Inf * ones(len(buses))
 236         Val = -Vau
 237         Vau[refs] = Va[refs]
 238         Val[refs] = Va[refs]
 239
 240         return Variable("Va", len(buses), Va, Val, Vau)
 241
 242
 243     def _get_voltage_magnitude_var(self, buses, generators):
 244         """ Returns the voltage magnitude variable set.
 245         """
 246         Vm = array([b.v_magnitude for b in buses])
 247
 248         # For buses with generators initialise Vm from gen data.
 249         for g in generators:
 250             Vm[g.bus._i] = g.v_magnitude
 251
 252         Vmin = array([b.v_min for b in buses])
 253         Vmax = array([b.v_max for b in buses])
 254
 255         return Variable("Vm", len(buses), Vm, Vmin, Vmax)
 256
 257
 258     def _get_pgen_var(self, generators, base_mva):
 259         """ Returns the generator active power set-point variable.
 260         """
 261         Pg = array([g.p / base_mva for g in generators])
 262
 263         Pmin = array([g.p_min / base_mva for g in generators])
 264         Pmax = array([g.p_max / base_mva for g in generators])
 265
 266         return Variable("Pg", len(generators), Pg, Pmin, Pmax)
 267
 268
 269     def _get_qgen_var(self, generators, base_mva):
 270         """ Returns the generator reactive power variable set.
 271         """
 272         Qg = array([g.q / base_mva for g in generators])
 273
 274         Qmin = array([g.q_min / base_mva for g in generators])
 275         Qmax = array([g.q_max / base_mva for g in generators])
 276
 277         return Variable("Qg", len(generators), Qg, Qmin, Qmax)
 278
 279     #--------------------------------------------------------------------------
 280     #  Constraints:
 281     #--------------------------------------------------------------------------
 282
 283     def _nln_constraints(self, nb, nl):
 284         """ Returns non-linear constraints for OPF.
 285         """
 286         Pmis = NonLinearConstraint("Pmis", nb)
 287         Qmis = NonLinearConstraint("Qmis", nb)
 288         Sf = NonLinearConstraint("Sf", nl)
 289         St = NonLinearConstraint("St", nl)
 290
 291         return Pmis, Qmis, Sf, St
 292
 293
 294     def _power_mismatch_dc(self, buses, generators, B, Pbusinj, base_mva):
 295         """ Returns the power mismatch constraint (B*Va + Pg = Pd).
 296         """
 297         nb, ng = len(buses), len(generators)
 298         # Negative bus-generator incidence matrix.
 299         gen_bus = array([g.bus._i for g in generators])
 300         neg_Cg = csr_matrix((-ones(ng), (gen_bus, range(ng))), (nb, ng))
 301
 302         Amis = hstack([B, neg_Cg], format="csr")
 303
 304         Pd = array([bus.p_demand for bus in buses])
 305         Gs = array([bus.g_shunt for bus in buses])
 306
 307         bmis = -(Pd - Gs) / base_mva - Pbusinj
 308
 309         return LinearConstraint("Pmis", Amis, bmis, bmis, ["Va", "Pg"])
 310
 311
 312     def _branch_flow_dc(self, branches, Bf, Pfinj, base_mva):
 313         """ Returns the branch flow limit constraint.  The real power flows
 314         at the from end the lines are related to the bus voltage angles by
 315         Pf = Bf * Va + Pfinj.
 316         """
 317         # Indexes of constrained lines.
 318         il = array([i for i,l in enumerate(branches) if 0.0 < l.rate_a < 1e10])
 319         lpf = -Inf * ones(len(il))
 320         rate_a = array([l.rate_a / base_mva for l in branches])
 321         upf = rate_a[il] - Pfinj[il]
 322         upt = rate_a[il] + Pfinj[il]
 323
 324         Pf = LinearConstraint("Pf",  Bf[il, :], lpf, upf, ["Va"])
 325         Pt = LinearConstraint("Pt", -Bf[il, :], lpf, upt, ["Va"])
 326
 327         return Pf, Pt
 328
 329
 330     def _const_pf_constraints(self, gn, base_mva):
 331         """ Returns a linear constraint enforcing constant power factor for
 332         dispatchable loads.
 333
 334         The power factor is derived from the original value of Pmin and either
 335         Qmin (for inductive loads) or Qmax (for capacitive loads). If both Qmin
 336         and Qmax are zero, this implies a unity power factor without the need
 337         for an additional constraint.
 338         """
 339         ivl = array([i for i, g in enumerate(gn)
 340                      if g.is_load and (g.q_min != 0.0 or g.q_max != 0.0)])
 341         vl = [gn[i] for i in ivl]
 342         nvl = len(vl)
 343
 344         ng = len(gn)
 345         Pg = array([g.p for g in vl]) / base_mva
 346         Qg = array([g.q for g in vl]) / base_mva
 347         Pmin = array([g.p_min for g in vl]) / base_mva
 348         Qmin = array([g.q_min for g in vl]) / base_mva
 349         Qmax = array([g.q_max for g in vl]) / base_mva
 350
 351         # At least one of the Q limits must be zero (corresponding to Pmax==0).
 352         for g in vl:
 353             if g.qmin != 0.0 and g.q_max != 0.0:
 354                 logger.error("Either Qmin or Qmax must be equal to zero for "
 355                 "each dispatchable load.")
 356
 357         # Initial values of PG and QG must be consistent with specified power
 358         # factor. This is to prevent a user from unknowingly using a case file
 359         # which would have defined a different power factor constraint under a
 360         # previous version which used PG and QG to define the power factor.
 361         Qlim = (Qmin == 0.0) * Qmax + (Qmax == 0.0) * Qmin
 362         if any( abs(Qg - Pg * Qlim / Pmin) > 1e-6 ):
 363             logger.error("For a dispatchable load, PG and QG must be "
 364                          "consistent with the power factor defined by "
 365                          "PMIN and the Q limits.")
 366
 367         # Make Avl, lvl, uvl, for lvl <= Avl * r_[Pg, Qg] <= uvl
 368         if nvl > 0:
 369             xx = Pmin
 370             yy = Qlim
 371             pftheta = arctan2(yy, xx)
 372             pc = sin(pftheta)
 373             qc = -cos(pftheta)
 374             ii = array([range(nvl), range(nvl)])
 375             jj = r_[ivl, ivl + ng]
 376             Avl = csr_matrix(r_[pc, qc], (ii, jj), (nvl, 2 * ng))
 377             lvl = zeros(nvl)
 378             uvl = lvl
 379         else:
 380             Avl = zeros((0, 2 * ng))
 381             lvl = array([])
 382             uvl = array([])
 383
 384         return LinearConstraint("vl", Avl, lvl, uvl, ["Pg", "Qg"])
 385
 386
 387     def _voltage_angle_diff_limit(self, buses, branches):
 388         """ Returns the constraint on the branch voltage angle differences.
 389         """
 390         nb = len(buses)
 391
 392         if not self.ignore_ang_lim:
 393             iang = [i for i, b in enumerate(branches)
 394                     if (b.ang_min and (b.ang_min > -360.0))
 395                     or (b.ang_max and (b.ang_max < 360.0))]
 396             iangl = array([i for i, b in enumerate(branches)
 397                      if b.ang_min is not None])[iang]
 398             iangh = array([i for i, b in enumerate(branches)
 399                            if b.ang_max is not None])[iang]
 400             nang = len(iang)
 401
 402             if nang > 0:
 403                 ii = range(nang) + range(nang)
 404                 jjf = array([b.from_bus._i for b in branches])[iang]
 405                 jjt = array([b.to_bus._i for b in branches])[iang]
 406                 jj = r_[jjf, jjt]
 407                 Aang = csr_matrix(r_[ones(nang), -ones(nang)], (ii, jj))
 408                 uang = Inf * ones(nang)
 409                 lang = -uang
 410                 lang[iangl] = array([b.ang_min * (pi / 180.0)
 411                                     for b in branches])[iangl]
 412                 uang[iangh] = array([b.ang_max * (pi / 180.0)
 413                                     for b in branches])[iangh]
 414             else:
 415 #                Aang = csr_matrix((0, nb), dtype=float64)
 416 #                lang = array([], dtype=float64)
 417 #                uang = array([], dtype=float64)
 418                 Aang = zeros((0, nb))
 419                 lang = array([])
 420                 uang = array([])
 421         else:
 422 #            Aang = csr_matrix((0, nb), dtype=float64)
 423 #            lang = array([], dtype=float64)
 424 #            uang = array([], dtype=float64)
 425 #            iang = array([], dtype=float64)
 426             Aang = zeros((0, nb))
 427             lang = array([])
 428             uang = array([])
 429
 430         return LinearConstraint("ang", Aang, lang, uang, ["Va"])
 431
 432
 433     def _pwl_gen_costs(self, generators, base_mva):
 434         """ Returns the basin constraints for piece-wise linear gen cost
 435         variables.  CCV cost formulation expressed as Ay * x <= by.
 436
 437         Based on makeAy.m from MATPOWER by C. E. Murillo-Sanchez, developed at
 438         PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
 439         information.
 440         """
 441         ng = len(generators)
 442         gpwl = [g for g in generators if g.pcost_model == PW_LINEAR]
 443 #        nq = len([g for g in gpwl if g.qcost_model is not None])
 444
 445         if self.dc:
 446             pgbas = 0        # starting index within x for active sources
 447             nq = 0           # number of Qg vars
 448 #            qgbas = None     # index of 1st Qg column in Ay
 449             ybas = ng        # starting index within x for y variables
 450         else:
 451             pgbas = 0
 452             nq = ng
 453 #            qgbas = ng + 1 # index of 1st Qg column in Ay
 454             ybas = ng + nq
 455
 456         # Number of extra y variables.
 457         ny = len(gpwl)
 458
 459         if ny == 0:
 460             return None, None
 461
 462         # Total number of cost points.
 463         nc = len([co for gn in gpwl for co in gn.p_cost])
 464 #        Ay = lil_matrix((nc - ny, ybas + ny))
 465         # Fill rows and then transpose.
 466         Ay = lil_matrix((ybas + ny, nc - ny))
 467         by = array([])
 468
 469         j = 0
 470         k = 0
 471         for i, g in enumerate(gpwl):
 472             # Number of cost points: segments = ns-1
 473             ns = len(g.p_cost)
 474
 475             p = array([x / base_mva for x, c in g.p_cost])
 476             c = array([c for x, c in g.p_cost])
 477             m = diff(c) / diff(p)        # Slopes for Pg (or Qg).
 478
 479             if 0.0 in diff(p):
 480                 raise ValueError, "Bad Pcost data: %s (%s)" % (p, g.name)
 481                 logger.error("Bad Pcost data: %s" % p)
 482
 483             b = m * p[:ns-1] - c[:ns-1] # rhs
 484             by = r_[by, b.T]
 485
 486 #            if i > ng:
 487 #                sidx = qgbas + (i-ng) - 1       # this was for a q cost
 488 #            else:
 489 #                sidx = pgbas + i - 1            # this was for a p cost
 490
 491             Ay[pgbas + i, k:k + ns - 1] = m
 492
 493             # FIXME: Repeat for Q costs.
 494
 495             # Now fill the y rows with -1's
 496             Ay[ybas + j, k:k + ns - 1] = -ones(ns-1)
 497
 498             k += (ns - 1)
 499             j += 1
 500
 501         y = Variable("y", ny)
 502
 503         # Transpose Ay since lil_matrix stores in rows.
 504         if self.dc:
 505             ycon = LinearConstraint("ycon", Ay.T, None, by, ["Pg", "y"])
 506         else:
 507             ycon = LinearConstraint("ycon", Ay.T, None, by, ["Pg", "Qg","y"])
 508
 509         return y, ycon
 510
 511 #------------------------------------------------------------------------------
 512 #  "UDOPF" class:
 513 #------------------------------------------------------------------------------
 514
 515 class UDOPF(OPF):
 516     """ The standard OPF formulation has no mechanism for completely shutting
 517     down generators which are very expensive to operate. Instead they are
 518     simply dispatched at their minimum generation limits. PYLON includes the
 519     capability to run an optimal power flow combined with a unit decommitment
 520     for a single time period, which allows it to shut down these expensive
 521     units and find a least cost commitment and dispatch.
 522
 523     Solves a combined unit decommitment and optimal power flow for a
 524     single time period. Uses an algorithm similar to dynamic programming.
 525     It proceeds through a sequence of stages, where stage N has N generators
 526     shut down, starting with N=0. In each stage, it forms a list of candidates
 527     (gens at their Pmin limits) and computes the cost with each one of them
 528     shut down. It selects the least cost case as the starting point for the
 529     next stage, continuing until there are no more candidates to be shut down
 530     or no more improvement can be gained by shutting something down.
 531
 532     Based on uopf.m from MATPOWER by Ray Zimmerman, developed at PSERC
 533     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
 534
 535     See also:
 536       - Ray Zimmerman, "MATPOWER User's Manual", MATPOWER, PSERC Cornell,
 537         version 3.2, U{http://www.pserc.cornell.edu/matpower/}, Sept, 2007
 538     """
 539
 540     def solve(self, solver_klass=None):
 541         """ Solves the combined unit decommitment / optimal power flow problem.
 542         """
 543         case = self.case
 544         generators = case.online_generators
 545
 546         logger.info("Solving OPF with unit de-commitment [%s]." % case.name)
 547
 548         t0 = time()
 549
 550         # 1. Begin at stage zero (N = 0), assuming all generators are on-line
 551         # with all limits in place. At most one generator shutdown per stage.
 552         i_stage = 0
 553
 554         # Check for sum(p_min) > total load, decommit as necessary.
 555         online = [g for g in generators if not g.is_load]
 556         online_vload = [g for g in generators if g.is_load]
 557
 558         # Total dispatchable load capacity.
 559         vload_capacity = sum([g.p_min for g in online_vload])
 560         # Total load capacity.
 561         load_capacity = sum([b.p_demand for b in case.buses]) - vload_capacity
 562
 563         # Minimum total online generation capacity.
 564         p_min_tot = sum([g.p_min for g in online])
 565
 566         # Shutdown the most expensive units until the minimum generation
 567         # capacity is less than the total load capacity.
 568         while p_min_tot > load_capacity:
 569             i_stage += 1
 570             logger.debug("De-commitment stage %d." % i_stage)
 571
 572             # Find generator with the maximum average cost at Pmin.
 573             avg_pmin_cost = [g.total_cost(g.p_min) / g.p_min for g in online]
 574             # Select at random from maximal generators with equal cost.
 575             g_idx, _ = fair_max(avg_pmin_cost)
 576             generator = online[g_idx]
 577
 578             logger.info("Shutting down generator [%s] to satisfy all "
 579                         "p_min limits." % generator.name)
 580
 581             # Shut down most expensive unit.
 582             generator.online = False
 583
 584             # Update minimum generation capacity for while loop.
 585             online = [g for g in case.online_generators if not g.is_load]
 586             p_min_tot = sum([g.p_min for g in online])
 587
 588         # 2. Solve a normal OPF and save the solution as the current best.
 589         solution = super(UDOPF, self).solve(solver_klass)
 590
 591         logger.debug("Initial system cost: $%.3f" % solution["f"])
 592
 593         if not solution["converged"] == True:
 594             logger.error("Non-convergent UDOPF [%s]." %
 595                          solution["output"]["message"])
 596             return solution
 597
 598         # 3. Go to the next stage, N = N + 1. Using the best solution from the
 599         # previous stage as the base case for this stage, ...
 600
 601         # Best case so far. A list of the on-line status of all generators.
 602         overall_online = [g.online for g in case.generators]
 603         # The objective function value is the total system cost.
 604         overall_cost = solution["f"]
 605
 606         # Best case for this stage.
 607         stage_online = overall_online
 608         stage_cost = overall_cost
 609
 610         # Shutdown at most one generator per stage.
 611         while True:
 612             # 4. Form a candidate list of generators with minimum
 613             # generation limits binding.
 614
 615             # Activate generators according to the stage best.
 616             for i, generator in enumerate(case.generators):
 617                 generator.online = stage_online[i]
 618
 619             # Get candidates for shutdown. Lagrangian multipliers are often
 620             # very small so we round to four decimal places.
 621             candidates = [g for g in case.online_generators if \
 622                           (round(g.mu_pmin, 4) > 0.0) and (g.p_min > 0.0)]
 623
 624             if len(candidates) == 0:
 625                 break
 626
 627             # Assume no improvement during this stage.
 628             done = True
 629
 630             i_stage += 1
 631             logger.debug("De-commitment stage %d." % i_stage)
 632
 633             for candidate in candidates:
 634                 # 5. For each generator on the candidate list, solve an OPF to
 635                 # find the total system cost with the generator shut down.
 636
 637                 # Activate generators according to the stage best.
 638                 for i, generator in enumerate(case.generators):
 639                     generator.online = stage_online[i]
 640
 641                 # Shutdown candidate generator.
 642                 candidate.online = False
 643
 644                 logger.debug("Solving OPF with generator '%s' shutdown." %
 645                     candidate.name)
 646
 647                 # Run OPF.
 648                 solution = super(UDOPF, self).solve(solver_klass)
 649
 650                 # Compare total system costs for improvement.
 651                 if solution["converged"] == True \
 652                     and (solution["f"] < overall_cost):
 653                     logger.debug("System cost improvement: $%.3f ($%.3f)" %
 654                                  (stage_cost - solution["f"], solution["f"]))
 655                     # 6. Replace the current best solution with this one if
 656                     # it has a lower cost.
 657                     overall_online = [g.online for g in case.generators]
 658                     overall_cost = solution["f"]
 659                     best_candidate = candidate
 660                     # Check for further decommitment.
 661                     done = False
 662                 else:
 663                     logger.debug("Candidate OPF failed [%s]." %
 664                                  solution["output"]["message"])
 665
 666                 # Reactivate the candidate before deactivating the next.
 667 #                candidate.online = True
 668
 669             if done:
 670                 # Decommits at this stage did not help.
 671                 break
 672             else:
 673                 # 7. If any of the candidate solutions produced an improvement,
 674                 # return to step 3.
 675
 676                 # Shutting something else down helps, so let's keep going.
 677                 logger.info("Shutting down generator '%s'.",
 678                             best_candidate.name)
 679
 680                 stage_online = overall_online
 681                 stage_cost = overall_cost
 682
 683         # 8. Use the best overall solution as the final solution.
 684         for i, generator in enumerate(case.generators):
 685             generator.online = overall_online[i]
 686
 687         # One final solve using the best case to ensure all results are
 688         # up-to-date.
 689         solution = super(UDOPF, self).solve(solver_klass)
 690
 691         logger.debug("UDOPF system cost: $%.3f" % solution["f"])
 692
 693         # Compute elapsed time and log it.
 694         elapsed = time() - t0
 695
 696         plural = "" if i_stage == 1 else "s"
 697         logger.info("Unit decommitment OPF solved in %.3fs (%d decommitment "
 698                     "stage%s)." % (elapsed, i_stage, plural))
 699
 700         return solution
 701
 702 #------------------------------------------------------------------------------
 703 #  "OPFModel" class:
 704 #------------------------------------------------------------------------------
 705
 706 class OPFModel(object):
 707     """ Defines a model for optimal power flow.
 708
 709     Based on @opf_model in MATPOWER by Ray Zimmerman, developed at PSERC
 710     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
 711     """
 712
 713     def __init__(self, case):
 714         #: Case to which the model relates.
 715         self.case = case
 716         #: Optimisation variables.
 717         self.vars = []
 718         #: Linear constraints.
 719         self.lin_constraints = []
 720         #: Non-linear constraints.
 721         self.nln_constraints = []
 722         #: User defined costs.
 723         self.costs = []
 724
 725
 726     @property
 727     def var_N(self):
 728         return sum([v.N for v in self.vars])
 729
 730
 731     def add_var(self, var):
 732         """ Adds a variable to the model.
 733         """
 734         if var.name in [v.name for v in self.vars]:
 735             logger.error("Variable set named '%s' already exists." % var.name)
 736             return
 737
 738         var.i1 = self.var_N
 739         var.iN = self.var_N + var.N - 1
 740         self.vars.append(var)
 741
 742
 743     def add_vars(self, vars):
 744         """ Adds a set of variables to the model.
 745         """
 746         for var in vars:
 747             self.add_var(var)
 748
 749
 750     def get_var(self, name):
 751         """ Returns the variable set with the given name.
 752         """
 753         for var in self.vars:
 754             if var.name == name:
 755                 return var
 756         else:
 757             raise ValueError
 758
 759
 760
 761     def get_var_N(self, name):
 762         """ Return the number of variables in the named set.
 763         """
 764         return self.get_var(name).N
 765
 766
 767     @property
 768     def nln_N(self):
 769         return sum([c.N for c in self.nln_constraints])
 770
 771
 772     @property
 773     def lin_N(self):
 774         return sum([c.N for c in self.lin_constraints])
 775
 776
 777     @property
 778     def lin_NS(self):
 779         return len(self.lin_constraints)
 780
 781
 782     def linear_constraints(self):
 783         """ Returns the linear constraints.
 784         """
 785         if self.lin_N == 0:
 786             return None, array([]), array([])
 787
 788         A = lil_matrix((self.lin_N, self.var_N), dtype=float64)
 789         l = -Inf * ones(self.lin_N)
 790         u = -l
 791
 792         for lin in self.lin_constraints:
 793             if lin.N:                   # non-zero number of rows to add
 794                 Ak = lin.A              # A for kth linear constrain set
 795                 i1 = lin.i1             # starting row index
 796                 iN = lin.iN             # ending row index
 797                 vsl = lin.vs            # var set list
 798                 kN = -1                 # initialize last col of Ak used
 799                 Ai = lil_matrix((lin.N, self.var_N), dtype=float64)
 800                 for v in vsl:
 801                     var = self.get_var(v)
 802                     j1 = var.i1         # starting column in A
 803                     jN = var.iN         # ending column in A
 804                     k1 = kN + 1         # starting column in Ak
 805                     kN = kN + var.N     # ending column in Ak
 806
 807                     if j1 == jN:
 808                         # FIXME: Single column slicing broken in lil.
 809                         for i in range(Ai.shape[0]):
 810                             Ai[i, j1] = Ak[i, k1]
 811                     else:
 812                         Ai[:, j1:jN + 1] = Ak[:, k1:kN + 1]
 813
 814                 A[i1:iN + 1, :] = Ai
 815                 l[i1:iN + 1] = lin.l
 816                 u[i1:iN + 1] = lin.u
 817
 818         return A.tocsr(), l, u
 819
 820
 821     def add_constraint(self, con):
 822         """ Adds a constraint to the model.
 823         """
 824         if isinstance(con, LinearConstraint):
 825             N, M = con.A.shape
 826             if con.name in [c.name for c in self.lin_constraints]:
 827                 logger.error("Constraint set named '%s' already exists."
 828                              % con.name)
 829                 return False
 830             else:
 831                 con.i1 = self.lin_N# + 1
 832                 con.iN = self.lin_N + N - 1
 833
 834                 nv = 0
 835                 for vs in con.vs:
 836                     nv = nv + self.get_var_N(vs)
 837                 if M != nv:
 838                     logger.error("Number of columns of A does not match number"
 839                         " of variables, A is %d x %d, nv = %d", N, M, nv)
 840                 self.lin_constraints.append(con)
 841         elif isinstance(con, NonLinearConstraint):
 842             N = con.N
 843             if con.name in [c.name for c in self.nln_constraints]:
 844                 logger.error("Constraint set named '%s' already exists."
 845                              % con.name)
 846                 return False
 847             else:
 848                 con.i1 = self.nln_N# + 1
 849                 con.iN = self.nln_N + N
 850                 self.nln_constraints.append(con)
 851         else:
 852             raise ValueError
 853
 854         return True
 855
 856
 857     def add_constraints(self, constraints):
 858         """ Adds constraints to the model.
 859         """
 860         for con in constraints:
 861             self.add_constraint(con)
 862
 863
 864     def get_lin_constraint(self, name):
 865         """ Returns the constraint set with the given name.
 866         """
 867         for c in self.lin_constraints:
 868             if c.name == name:
 869                 return c
 870         else:
 871             raise ValueError
 872
 873
 874     def get_nln_constraint(self, name):
 875         """ Returns the constraint set with the given name.
 876         """
 877         for c in self.nln_constraints:
 878             if c.name == name:
 879                 return c
 880         else:
 881             raise ValueError
 882
 883     @property
 884     def cost_N(self):
 885         return sum([c.N for c in self.costs])
 886
 887
 888     def get_cost_params(self):
 889         """ Returns the cost parameters.
 890         """
 891         return [c.params for c in self.costs]
 892
 893 #------------------------------------------------------------------------------
 894 #  "_Set" class:
 895 #------------------------------------------------------------------------------
 896
 897 class _Set(_Named):
 898
 899     def __init__(self, name, N):
 900         #: String identifier.
 901         self.name = name
 902
 903         #: Starting index.
 904         self.i0 = 0
 905
 906         #: Ending index.
 907         self.iN = 0
 908
 909         #: Number in set.
 910         self.N = N
 911
 912         #: Number of sets.
 913         self.NS = 0
 914
 915         #: Ordered list of sets.
 916         self.order = []
 917
 918 #------------------------------------------------------------------------------
 919 #  "Variable" class:
 920 #------------------------------------------------------------------------------
 921
 922 class Variable(_Set):
 923     """ Defines a set of variables.
 924     """
 925
 926     def __init__(self, name, N, v0=None, vl=None , vu=None):
 927         """ Initialises a new Variable instance.
 928         """
 929         super(Variable, self).__init__(name, N)
 930
 931         #: Initial value of the variables. Zero by default.
 932         if v0 is None:
 933             self.v0 = zeros(N)
 934         else:
 935             self.v0 = v0
 936
 937         #: Lower bound on the variables. Unbounded below be default.
 938         if vl is None:
 939             self.vl = -Inf * ones(N)
 940         else:
 941             self.vl = vl
 942
 943         #: Upper bound on the variables. Unbounded above by default.
 944         if vu is None:
 945             self.vu = Inf * ones(N)
 946         else:
 947             self.vu = vu
 948
 949 #------------------------------------------------------------------------------
 950 #  "LinearConstraint" class:
 951 #------------------------------------------------------------------------------
 952
 953 class LinearConstraint(_Set):
 954     """ Defines a set of linear constraints.
 955     """
 956
 957     def __init__(self, name, AorN, l=None, u=None, vs=None):
 958         N, _ = AorN.shape
 959
 960         super(LinearConstraint, self).__init__(name, N)
 961
 962         self.A = AorN
 963         self.l = -Inf * ones(N) if l is None else l
 964         self.u =  Inf * ones(N) if u is None else u
 965
 966         # Varsets.
 967         self.vs = [] if vs is None else vs
 968
 969         if (self.l.shape[0] != N) or (self.u.shape[0] != N):
 970             logger.error("Sizes of A, l and u must match.")
 971
 972 #------------------------------------------------------------------------------
 973 #  "NonLinearConstraint" class:
 974 #------------------------------------------------------------------------------
 975
 976 class NonLinearConstraint(_Set):
 977     """ Defines a set of non-linear constraints.
 978     """
 979     pass
 980
 981 #------------------------------------------------------------------------------
 982 #  "Cost" class:
 983 #------------------------------------------------------------------------------
 984
 985 class Cost(_Set):
 986     def __init__(self):
 987         self.N = None
 988         self.H = None
 989         self.Cw = None
 990         self.dd = None
 991         self.rh = None
 992         self.kk = None
 993         self.mm = None
 994         self.vs = None
 995         self.params = None
 996
 997 # EOF -------------------------------------------------------------------------