Computes 2nd derivatives of power injection w.r.t. voltage.
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
88151
def PartialDynamicSystem(self, ieq, variable):
        
        if ieq == 0:
            # U2-U1=signal
            if variable == self.physical_nodes[0].variable:
                # U1 is output
                # U1=U2-signal
                return [WeightedSum([self.physical_nodes[1].variable, self.voltage_signal], variable, [1, -1])]
            elif variable == self.physical_nodes[1].variable:
                # U2 is output
                # U2=U1+signal
                return [WeightedSum([self.physical_nodes[0].variable, self.voltage_signal], variable, [1, 1])]
602883
def derivatives(self, x, y, n_sersic, R_sersic, k_eff, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        e = abs(1. - q)
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        x_, y_ = self._coord_transf(x, y, q, phi_G, center_x, center_y)
        f_x_prim, f_y_prim = self.sersic.derivatives(x_, y_, n_sersic, R_sersic, k_eff)
        f_x_prim *= np.sqrt(1 - e)
        f_y_prim *= np.sqrt(1 + e)
        f_x = cos_phi*f_x_prim-sin_phi*f_y_prim
        f_y = sin_phi*f_x_prim+cos_phi*f_y_prim
        return f_x, f_y
725069
public static double PowerSeries(double a, double b, double x) {
        double s, t, u, v, n, t1, z, ai;

        ai = 1.0 / a;
        u = (1.0 - b) * x;
        v = u / (a + 1.0);
        t1 = v;
        t = u;
        n = 2.0;
        s = 0.0;
        z = Constants.DoubleEpsilon * ai;
        while (Math.abs(v) > z) {
            u = (n - b) * x / n;
            t *= u;
            v = t / (a + n);
            s += v;
            n += 1.0;
        }
        s += t1;
        s += ai;

        u = a * Math.log(x);
        if ((a + b) < Gamma.GammaMax && Math.abs(u) < Constants.LogMax) {
            t = Gamma.Function(a + b) / (Gamma.Function(a) * Gamma.Function(b));
            s = s * t * Math.pow(x, a);
        } else {
            t = Gamma.Log(a + b) - Gamma.Log(a) - Gamma.Log(b) + u + Math.log(s);
            if (t < Constants.LogMin) s = 0.0;
            else s = Math.exp(t);
        }
        return s;
    }
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
601528
def post_process_data(self, voltage=None, incremental_capacity=None,
                          voltage_step=None):
        """"""

        logging.debug("post-processing data")

        if voltage is None:
            voltage = self.voltage_processed
            incremental_capacity = self.incremental_capacity
            voltage_step = self.voltage_inverted_step

        if self.post_smoothing:
            logging.debug(" - post smoothing (gaussian)")
            logging.debug(f"    * using voltage fwhm: {self.voltage_fwhm}")
            points_fwhm = int(self.voltage_fwhm / voltage_step)

            sigma = np.amax([1, points_fwhm / 2])

            self.incremental_capacity = gaussian_filter1d(
                incremental_capacity, sigma=sigma, order=self.gaussian_order,
                mode=self.gaussian_mode,
                cval=self.gaussian_cval,
                truncate=self.gaussian_truncate
            )

        if self.normalize:
            logging.debug(" - normalizing")
            area = simps(incremental_capacity, voltage)
            self.incremental_capacity = incremental_capacity * self.normalizing_factor / abs(area)

        fixed_range = False
        if isinstance(self.fixed_voltage_range, np.ndarray):
            fixed_range = True
        else:
            if self.fixed_voltage_range:
                fixed_range = True
        if fixed_range:
            logging.debug(" - using fixed voltage range (interpolating)")
            v1, v2, number_of_points = self.fixed_voltage_range
            v = np.linspace(v1, v2, number_of_points)
            f = interp1d(x=self.voltage_processed, y=self.incremental_capacity,
                         kind=self.interpolation_method, bounds_error=False,
                         fill_value=np.NaN)

            self.incremental_capacity = f(v)
            self.voltage_processed = v
366856
def _deriv_logaddexp2(x1, x2):
    """"""
    y1 = np.exp2(x1)
    y2 = np.exp2(x2)
    df_dx1 = y1 / (y1 + y2)
    df_dx2 = y2 / (y1 + y2)
    return np.vstack([df_dx1, df_dx2]).T
80294
def sub_v2(v, w):
    

    if type(w) is float or type(w) is int:
        return Vec2(v.x - w, v.y - w)
    else:
        return Vec2(v.x - w.x, v.y - w.y)
744886
def dSbus_dV(Y, V):
    
    I = Y * V

    diagV = spdiag(V)
    diagIbus = spdiag(I)
    diagVnorm = spdiag(div(V, abs(V))) # Element-wise division.

    dS_dVm = diagV * conj(Y * diagVnorm) + conj(diagIbus) * diagVnorm
    dS_dVa = 1j * diagV * conj(diagIbus - Y * diagV)

    return dS_dVm, dS_dVa
349473
def d2AIbr_dV2(dIbr_dVa, dIbr_dVm, Ibr, Ybr, V, lam):
    
    diaglam = spdiag(lam)
    diagIbr_conj = spdiag(conj(Ibr))

    Iaa, Iav, Iva, Ivv = d2Ibr_dV2(Ybr, V, diagIbr_conj * lam)

    Haa = 2 * ( Iaa + dIbr_dVa.T * diaglam * conj(dIbr_dVa) ).real()
    Hva = 2 * ( Iva + dIbr_dVm.T * diaglam * conj(dIbr_dVa) ).real()
    Hav = 2 * ( Iav + dIbr_dVa.T * diaglam * conj(dIbr_dVm) ).real()
    Hvv = 2 * ( Ivv + dIbr_dVm.T * diaglam * conj(dIbr_dVm) ).real()

    return Haa, Hav, Hva, Hvv
1076524
def _hessfcn(self, x, lmbda):
        
        Pgen = x[self._Pg.i1:self._Pg.iN + 1] # Active generation in p.u.
        Qgen = x[self._Qg.i1:self._Qg.iN + 1] # Reactive generation in p.u.

        for i, g in enumerate(self._gn):
            g.p = Pgen[i] * self._base_mva # active generation in MW
            g.q = Qgen[i] * self._base_mva # reactive generation in MVAr

        Vang = x[self._Va.i1:self._Va.iN + 1]
        Vmag = x[self._Vm.i1:self._Vm.iN + 1]
        V = Vmag * exp(1j * Vang)
        nxtra = self._nxyz - 2 * self._nb

        #------------------------------------------------------------------
        #  Evaluate d2f.
        #------------------------------------------------------------------

        d2f = self._d2f(x) * self.opt["cost_mult"]
        # TODO: Generalised cost model.

        #------------------------------------------------------------------
        #  Evaluate Hessian of power balance constraints.
        #------------------------------------------------------------------

        nlam = len(lmbda["eqnonlin"]) / 2
        lamP = lmbda["eqnonlin"][:nlam]
        lamQ = lmbda["eqnonlin"][nlam:nlam + nlam]
        Gpaa, Gpav, Gpva, Gpvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamP)
        Gqaa, Gqav, Gqva, Gqvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamQ)

        d2G = vstack([
            hstack([
                vstack([hstack([Gpaa, Gpav]),
                        hstack([Gpva, Gpvv])]).real +
                vstack([hstack([Gqaa, Gqav]),
                        hstack([Gqva, Gqvv])]).imag,
                csr_matrix((2 * self._nb, nxtra))]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        #------------------------------------------------------------------
        #  Evaluate Hessian of flow constraints.
        #------------------------------------------------------------------

        nmu = len(lmbda["ineqnonlin"]) / 2
        muF = lmbda["ineqnonlin"][:nmu]
        muT = lmbda["ineqnonlin"][nmu:nmu + nmu]
        if self.flow_lim == "I":
            dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
                self.om.case.dIbr_dV(self._Yf, self._Yt, V)
            Hfaa, Hfav, Hfva, Hfvv = \
                self.om.case.d2AIbr_dV2(dIf_dVa, dIf_dVm, If, self._Yf, V, muF)
            Htaa, Htav, Htva, Htvv = \
                self.om.case.d2AIbr_dV2(dIt_dVa, dIt_dVm, It, self._Yt, V, muT)
        else:
            f = [e.from_bus._i for e in self._ln]
            t = [e.to_bus._i for e in self._ln]
            # Line-bus connection matrices.
            Cf = csr_matrix((ones(self._nl), (range(self._nl), f)), (self._nl, self._nb))
            Ct = csr_matrix((ones(self._nl), (range(self._nl), t)), (self._nl, self._nb))
            dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
                self.om.case.dSbr_dV(self._Yf, self._Yt, V)
            if self.flow_lim == PFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
                                            Sf.real(), Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
                                            St.real(), Ct, self._Yt, V, muT)
            elif self.flow_lim == SFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(
                        dSf_dVa, dSf_dVm, Sf, Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(
                        dSt_dVa, dSt_dVm, St, Ct, self._Yt, V, muT)
            else:
                raise ValueError

        d2H = vstack([
            hstack([
                vstack([hstack([Hfaa, Hfav]),
                        hstack([Hfva, Hfvv])]) +
                vstack([hstack([Htaa, Htav]),
                        hstack([Htva, Htvv])]),
                csr_matrix((2 * self._nb, nxtra))
            ]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        return d2f + d2G + d2H
1292017
def partial_derivative(self, X, y=0):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 0:
            return V

        else:
            A = np.power(V, -self.theta - 1)
            B = np.power(V, -self.theta) + np.power(U, -self.theta) - 1
            h = np.power(B, (-1 - self.theta) / self.theta)
            return np.multiply(A, h) - y
1063410
def _deriv_logaddexp(x1, x2):
    """"""
    y1 = np.exp(x1)
    y2 = np.exp(x2)
    df_dx1 = y1 / (y1 + y2)
    df_dx2 = y2 / (y1 + y2)
    return np.vstack([df_dx1, df_dx2]).T
88006
def calculate_dependent_values(network):
    """"""

    apply_line_types(network)
    apply_transformer_types(network)

    network.lines["v_nom"] = network.lines.bus0.map(network.buses.v_nom)

    network.lines["x_pu"] = network.lines.x/(network.lines.v_nom**2)
    network.lines["r_pu"] = network.lines.r/(network.lines.v_nom**2)
    network.lines["b_pu"] = network.lines.b*network.lines.v_nom**2
    network.lines["g_pu"] = network.lines.g*network.lines.v_nom**2
    network.lines["x_pu_eff"] = network.lines["x_pu"]
    network.lines["r_pu_eff"] = network.lines["r_pu"]


    #convert transformer impedances from base power s_nom to base = 1 MVA
    network.transformers["x_pu"] = network.transformers.x/network.transformers.s_nom
    network.transformers["r_pu"] = network.transformers.r/network.transformers.s_nom
    network.transformers["b_pu"] = network.transformers.b*network.transformers.s_nom
    network.transformers["g_pu"] = network.transformers.g*network.transformers.s_nom
    network.transformers["x_pu_eff"] = network.transformers["x_pu"]* network.transformers["tap_ratio"]
    network.transformers["r_pu_eff"] = network.transformers["r_pu"]* network.transformers["tap_ratio"]

    apply_transformer_t_model(network)

    network.shunt_impedances["v_nom"] = network.shunt_impedances["bus"].map(network.buses.v_nom)
    network.shunt_impedances["b_pu"] = network.shunt_impedances.b*network.shunt_impedances.v_nom**2
    network.shunt_impedances["g_pu"] = network.shunt_impedances.g*network.shunt_impedances.v_nom**2
337604
def derivatives(self, x, y, Rs, theta_Rs, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        e = min(abs(1. - q), 0.99)
        xt1 = (cos_phi*x_shift+sin_phi*y_shift)*np.sqrt(1 - e)
        xt2 = (-sin_phi*x_shift+cos_phi*y_shift)*np.sqrt(1 + e)
        R_ = np.sqrt(xt1**2 + xt2**2)
        rho0_input = self.nfw._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs)
        if Rs < 0.0000001:
            Rs = 0.0000001
        f_x_prim, f_y_prim = self.nfw.nfwAlpha(R_, Rs, rho0_input, xt1, xt2)
        f_x_prim *= np.sqrt(1 - e)
        f_y_prim *= np.sqrt(1 + e)
        f_x = cos_phi*f_x_prim-sin_phi*f_y_prim
        f_y = sin_phi*f_x_prim+cos_phi*f_y_prim
        return f_x, f_y
47940
protected double derivative(int i, NumberVector v) {
    final int dim = v.getDimensionality();
    if(dim == 1) {
      return 0.;
    }
    // Adjust for boundary conditions, as per the article:
    i = (i == 0) ? 1 : (i == dim - 1) ? dim - 2 : i;
    return (v.doubleValue(i) - v.doubleValue(i - 1) + (v.doubleValue(i + 1) - v.doubleValue(i - 1)) * .5) * .5;
  }
579980
def _q_iteration(self, Q, Bpp_solver, Vm, Va, pq):
        
        dVm = -Bpp_solver.solve(Q)

        # Update voltage.
        Vm[pq] = Vm[pq] + dVm
        V = Vm * exp(1j * Va)

        return V, Vm, Va