Fast Decoupled Power Flow execution: Implement g(y)
968969
def _compile_fdpf(self):
        
        string = '"""\n'
        string += 'system.dae.init_g()\n'
        for pflow, gcall, call in zip(self.pflow, self.gcall, self.gcalls):
            if pflow and gcall:
                string += call
        string += 'system.dae.reset_small_g()\n'
        string += '\n'
        string += '"""'
        self.fdpf = compile(eval(string), '', 'exec')
399167
def solveNetwork(self, Xgen, augYbus_solver, gbus):
        
        generators = self.dyn_generators
        j = 0 + 1j

        ng = len(gbus)
        Igen = zeros(ng)

        s = len(augYbus_solver)
        Ig = zeros(s)

        # Define generator types.
        typ1 = [g._i for g in generators if g.model == CLASSICAL]
        typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]

        # Generator type 1: classical model
        delta = Xgen[typ1, 0]
        Eq_tr = Xgen[typ1, 2]

        xd_tr = array([g.xd_tr for g in generators])[typ1]

        # Calculate generator currents
        Igen[typ1] = (Eq_tr * exp(j * delta)) / (j * xd_tr)

        # Generator type 2: 4th order model
        delta = Xgen[typ2, 0]
        Eq_tr = Xgen[typ2, 2]
        Ed_tr = Xgen[typ2, 3]

        xd_tr = array([g.xd_tr for g in generators])[typ2] # Pgen(type2,8)

        # Calculate generator currents. (Padiyar, p.417.)
        Igen[typ2] = (Eq_tr + j * Ed_tr) * exp(j * delta) / (j * xd_tr)

        # Calculations --------------------------------------------------------

        # Generator currents
        Ig[gbus] = Igen

        # Calculate network voltages: U = Y/Ig
        U = augYbus_solver.solve(Ig)

        return U
989127
def define_passive_branch_flows_with_kirchhoff(network,snapshots,skip_vars=False):
    

    for sub_network in network.sub_networks.obj:
        find_tree(sub_network)
        find_cycles(sub_network)

        #following is necessary to calculate angles post-facto
        find_bus_controls(sub_network)
        if len(sub_network.branches_i()) > 0:
            calculate_B_H(sub_network)

    passive_branches = network.passive_branches()

    if not skip_vars:
        network.model.passive_branch_p = Var(list(passive_branches.index), snapshots)

    cycle_index = []
    cycle_constraints = {}

    for subnetwork in network.sub_networks.obj:

        branches = subnetwork.branches()
        attribute = "r_pu_eff" if network.sub_networks.at[subnetwork.name,"carrier"] == "DC" else "x_pu_eff"

        sub_network_cycle_index, sub_network_cycle_constraints = define_sub_network_cycle_constraints( subnetwork,
                                                                                                       snapshots,
                                                                                                       network.model.passive_branch_p, attribute)

        cycle_index.extend( sub_network_cycle_index)
        cycle_constraints.update( sub_network_cycle_constraints)

    l_constraint(network.model, "cycle_constraints", cycle_constraints,
                 cycle_index, snapshots)
398726
private Set<Integer> getSuccessors(Graph g, Integer vertex) {
        if (g instanceof DirectedGraph) {
            DirectedGraph<?> dg = (DirectedGraph)g;
            return dg.successors(vertex);
        }
        else
            return Collections.<Integer>emptySet();
    }
966088
def d2Sbus_dV2(Ybus, V, lam):
    
    n = len(V)
    Ibus = Ybus * V
    diaglam = spdiag(lam)
    diagV = spdiag(V)

    A = spmatrix(mul(lam, V), range(n), range(n))
    B = Ybus * diagV
    C = A * conj(B)
    D = Ybus.H * diagV
    E = conj(diagV) * (D * diaglam - spmatrix(D*lam, range(n), range(n)))
    F = C - A * spmatrix(conj(Ibus), range(n), range(n))
    G = spmatrix(div(matrix(1.0, (n, 1)), abs(V)), range(n), range(n))

    Gaa = E + F
    Gva = 1j * G * (E - F)
    Gav = Gva.T
    Gvv = G * (C + C.T) * G

    return Gaa, Gav, Gva, Gvv
1233930
def run(self):
        
        ret = False
        system = self.system
        config = self.config
        dae = self.system.dae

        # maxit = config.maxit
        # tol = config.tol

        if system.pflow.solved is False:
            logger.warning('Power flow not solved. Simulation cannot continue.')
            return ret
        t0, _ = elapsed()
        t1 = t0

        self.streaming_init()

        logger.info('')
        logger.info('-> Time Domain Simulation: {} method, t={} s'
                    .format(self.config.method, self.config.tf))

        self.load_pert()

        self.run_step0()

        config.qrtstart = time()

        while self.t < config.tf:
            self.check_fixed_times()
            self.calc_time_step()

            if self.callpert is not None:
                self.callpert(self.t, self.system)

            if self.h == 0:
                break
            # progress time and set time in dae
            self.t += self.h
            dae.t = self.t

            # backup actual variables
            self.x0 = matrix(dae.x)
            self.y0 = matrix(dae.y)
            self.f0 = matrix(dae.f)

            # apply fixed_time interventions and perturbations
            self.event_actions()

            # reset flags used in each step
            self.err = 1
            self.niter = 0
            self.convergence = False

            self.implicit_step()

            if self.convergence is False:
                try:
                    self.restore_values()
                    continue
                except ValueError:
                    self.t = config.tf
                    ret = False
                    break

            self.step += 1
            self.compute_flows()
            system.varout.store(self.t, self.step)
            self.streaming_step()

            # plot variables and display iteration status
            perc = max(min((self.t - config.t0) / (config.tf - config.t0) * 100, 100), 0)

            # show iteration info every 30 seconds or every 20%

            t2, _ = elapsed(t1)
            if t2 - t1 >= 30:
                t1 = t2
                logger.info(' ({:.0f}%) time = {:.4f}s, step = {}, niter = {}'
                            .format(100 * self.t / config.tf, self.t, self.step,
                                    self.niter))

            if perc > self.next_pc or self.t == config.tf:
                self.next_pc += 20
                logger.info(' ({:.0f}%) time = {:.4f}s, step = {}, niter = {}'
                            .format(100 * self.t / config.tf, self.t, self.step, self.niter))

            # compute max rotor angle difference
            # diff_max = anglediff()

            # quasi-real-time check and wait
            rt_end = config.qrtstart + (self.t - config.t0) * config.kqrt

            if config.qrt:
                # the ending time has passed
                if time() - rt_end > 0:
                    # simulation is too slow
                    if time() - rt_end > config.kqrt:
                        logger.debug('Simulation over-run at t={:4.4g} s.'.format(self.t))
                # wait to finish
                else:
                    self.headroom += (rt_end - time())
                    while time() - rt_end < 0:
                        sleep(1e-5)

        if config.qrt:
            logger.debug('RT headroom time: {} s.'.format(str(self.headroom)))

        if self.t != config.tf:
            logger.error('Reached minimum time step. Convergence is not likely.')
            ret = False
        else:
            ret = True

        if system.config.dime_enable:
            system.streaming.finalize()

        _, s = elapsed(t0)

        if ret is True:
            logger.info(' Time domain simulation finished in {:s}.'.format(s))
        else:
            logger.info(' Time domain simulation failed in {:s}.'.format(s))

        self.success = ret
        self.dump_results(success=self.success)

        return ret
107697
public final void pshufd(XMMRegister dst, XMMRegister src, Immediate imm8)
  {
    emitX86(INST_PSHUFD, dst, src, imm8);
  }
1253992
public void mul(int i, int j, double v)
    {
        consumer.set(i, j, supplier.get(i, j) * v);
    }
208523
private double optimize(int k0, int kmax, double sumx, double sumx2, double xp, double yp, double sumxy, double sumy) {
    int k = kmax - k0 + 1;
    return (sumxy - xp * sumy - yp * sumx + k * xp * yp) / (sumx2 - 2 * sumx * xp + k * xp * xp);
    // return (-xp * yp * + yp * sumx - sumxy + xp * sumy) / (-xp * xp *
    // kmax - sumx2 + 2 * xp * sumx);
  }
1116194
def performAction(self, action):
        
        gs = [g for g in self.case.online_generators if g.bus.type !=REFERENCE]

        assert len(action) == len(gs)

        logger.info("Action: %s" % list(action))

        # Set the output of each (non-reference) generator.
        for i, g in enumerate(gs):
            g.p = action[i]

        # Compute power flows and slack generator set-point.
        NewtonPF(self.case, verbose=False).solve()
        #FastDecoupledPF(self.case, verbose=False).solve()

        # Store all generator set-points (only used for plotting).
        self._Pg[:, self._step] = [g.p for g in self.case.online_generators]

        # Apply the next load profile value to the original demand at each bus.
        if self._step != len(self.profile) - 1:
            pq_buses = [b for b in self.case.buses if b.type == PQ]
            for i, b in enumerate(pq_buses):
                b.p_demand = self._Pd0[i] * self.profile[self._step + 1]

        self._step += 1

        logger.info("Entering step %d." % self._step)
1129129
def par_y0bstep(i):
    
    """

    global mp_Y0
    mp_Y0[i] = 1/mp_rho*mp_S + mp_DX[i] + mp_U0[i] + mp_b
321489
def _compute_value(power, wg):
    """"""
    if power not in wg:
        p1, p2 = power
        # y power
        if p1 == 0:
            yy = wg[(0, -1)]
            wg[power] = numpy.power(yy, p2 / 2).sum() / len(yy)
        # x power
        else:
            xx = wg[(-1, 0)]
            wg[power] = numpy.power(xx, p1 / 2).sum() / len(xx)
    return wg[power]
900145
function alterneigh(rad, i, b, g, r) {

    var j;
    var k;
    var lo;
    var hi;
    var a;
    var m;
    var p;

    lo = i - rad;
    if (lo < -1) lo = -1;

    hi = i + rad;
    if (hi > netsize) hi = netsize;

    j = i + 1;
    k = i - 1;
    m = 1;

    while ((j < hi) || (k > lo)) {
      a = radpower[m++];

      if (j < hi) {
        p = network[j++];

        try {
          p[0] -= (a * (p[0] - b)) / alpharadbias;
          p[1] -= (a * (p[1] - g)) / alpharadbias;
          p[2] -= (a * (p[2] - r)) / alpharadbias;
        } catch (e) { } // prevents 1.3 miscompilation
      }

      if (k > lo) {
        p = network[k--];

        try {
          p[0] -= (a * (p[0] - b)) / alpharadbias;
          p[1] -= (a * (p[1] - g)) / alpharadbias;
          p[2] -= (a * (p[2] - r)) / alpharadbias;
        } catch (e) { }
      }
    }
  }
89385
def _dinic_step(graph, capacity, lev, flow, u, target, limit):
    
    if limit <= 0:
        return 0
    if u == target:
        return limit
    val = 0
    for v in graph[u]:
        residual = capacity[u][v] - flow[u][v]
        if lev[v] == lev[u] + 1 and residual > 0:
            z = min(limit, residual)
            aug = _dinic_step(graph, capacity, lev, flow, v, target, z)
            flow[u][v] += aug
            flow[v][u] -= aug
            val += aug
            limit -= aug
    if val == 0:
        lev[u] = None         # remove unreachable node
    return val
32749
function (dx, dy) {

    	    var phynode_id = this.id;

	    // Move box positions.
    	    var assoc_phynodes = get_descendant_phynodes(phynode_id);
    	    for( var si = 0; si < assoc_phynodes.length; si++ ){
		var mshp = assoc_phynodes[si];
		mshp.move_y(dy);
		//ll('mshp['+si+']:'+' oy: '+mshp.start_y+', dy:'+dy);
    	    }

	    // Collect subtree edges for next bit.
	    var dimmable_subtree = {};
	    var subtree_edges = get_descendant_connections(phynode_id);
	    for( var se = 0; se < subtree_edges.length; se++ ){
		var ste = subtree_edges[se];
		dimmable_subtree[ste.id] = true;
	    }

	    // Update connections; keep subtree dimmed while in transit.
            for (var i = connections.length; i--;) {
		var conn = connections[i];
		if( dimmable_subtree[conn.id] ){
		    conn.update('dim');		
		}else{
		    conn.update();		
		}
            }
            paper.safari();
	}
1287251
def _run_opf(self, P, q, AA, ll, uu, xmin, xmax, x0, opt):
        
        AAcvx = tocvx(AA)
        nx = x0.shape[0] # number of variables
        # add var limits to linear constraints
        eyex = spmatrix(1.0, range(nx), range(nx))
        A = sparse([eyex, AAcvx])
        l = matrix([-xmin, ll])
        u = matrix([ xmax, uu])

        Ae, be, Ai, bi = split_linear_constraints(A, l, u)


        if len(P.V) > 0:
            cvx_sol = solvers.qp(P, q, Ai, bi, Ae, be, self.solver, {"x": x0})
        else:
            cvx_sol = solvers.lp(q, Ai, bi, Ae, be, self.solver, {"x": x0})

        return cvx_sol
175649
def _power_mismatch_dc(self, buses, generators, B, Pbusinj, base_mva):
        
        nb, ng = len(buses), len(generators)
        # Negative bus-generator incidence matrix.
        gen_bus = array([g.bus._i for g in generators])
        neg_Cg = csr_matrix((-ones(ng), (gen_bus, range(ng))), (nb, ng))

        Amis = hstack([B, neg_Cg], format="csr")

        Pd = array([bus.p_demand for bus in buses])
        Gs = array([bus.g_shunt for bus in buses])

        bmis = -(Pd - Gs) / base_mva - Pbusinj

        return LinearConstraint("Pmis", Amis, bmis, bmis, ["Va", "Pg"])