Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
dae.y[self.isq], self.xmu)) - mul(dae.x[self.ird], self.rr)
dae.g[self.vrq] = -dae.y[self.vrq] - mul(
dae.x[self.irq], self.rr) - mul(
1 - dae.x[self.omega_m],
mul(dae.x[self.ird], self.x1) + mul(dae.y[self.isd], self.xmu))
dae.g[self.vsd] = -dae.y[self.vsd] - mul(dae.y[self.v],
sin(dae.y[self.a]))
dae.g[self.vsq] = -dae.y[self.vsq] + mul(dae.y[self.v],
cos(dae.y[self.a]))
dae.g[self.vref] = self.vref0 - dae.y[self.vref]
dae.g[self.pwa] = mmax(mmin(2 * dae.x[self.omega_m] - 1, 1),
0) - dae.y[self.pwa]
dae.hard_limit(self.pwa, 0, 1)
dae.g[self.pw] = -dae.y[self.pw] + mul(
0.5, dae.y[self.cp], self.ngen, pi, self.rho, (self.R)**2,
(self.Vwn)**3, div(1, self.mva_mega), (dae.x[self.vw])**3)
dae.g[self.cp] = -dae.y[self.cp] + mul(
-1.1 + mul(25.52, div(1, dae.y[self.ilamb])) + mul(
-0.088, dae.x[self.theta_p]),
exp(mul(-12.5, div(1, dae.y[self.ilamb]))))
dae.g[self.lamb] = -dae.y[self.lamb] + mul(
4, self.R, self.fn, self.ngb, dae.x[self.omega_m], pi,
div(1, self.Vwn), div(1, self.npole), div(1, dae.x[self.vw]))
dae.g[self.ilamb] = div(
1,
div(1, dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p])) + mul(
-0.035, div(1, 1 +
(dae.x[self.theta_p])**3))) - dae.y[self.ilamb]
dae.g += spmatrix(
mul(
self.u0, -mul(dae.x[self.ird], dae.y[self.vrd]) - mul(
dae.g[self.vsq] = -dae.y[self.vsq] + mul(dae.y[self.v],
cos(dae.y[self.a]))
dae.g[self.vref] = self.vref0 - dae.y[self.vref]
dae.g[self.pwa] = mmax(mmin(2 * dae.x[self.omega_m] - 1, 1),
0) - dae.y[self.pwa]
dae.hard_limit(self.pwa, 0, 1)
dae.g[self.pw] = -dae.y[self.pw] + mul(
0.5, dae.y[self.cp], self.ngen, pi, self.rho, (self.R)**2,
(self.Vwn)**3, div(1, self.mva_mega), (dae.x[self.vw])**3)
dae.g[self.cp] = -dae.y[self.cp] + mul(
-1.1 + mul(25.52, div(1, dae.y[self.ilamb])) + mul(
-0.088, dae.x[self.theta_p]),
exp(mul(-12.5, div(1, dae.y[self.ilamb]))))
dae.g[self.lamb] = -dae.y[self.lamb] + mul(
4, self.R, self.fn, self.ngb, dae.x[self.omega_m], pi,
div(1, self.Vwn), div(1, self.npole), div(1, dae.x[self.vw]))
dae.g[self.ilamb] = div(
1,
div(1, dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p])) + mul(
-0.035, div(1, 1 +
(dae.x[self.theta_p])**3))) - dae.y[self.ilamb]
dae.g += spmatrix(
mul(
self.u0, -mul(dae.x[self.ird], dae.y[self.vrd]) - mul(
dae.x[self.irq], dae.y[self.vrq]) - mul(
dae.y[self.isd], dae.y[self.vsd]) - mul(
dae.y[self.isq], dae.y[self.vsq])), self.a,
[0] * self.n, (dae.m, 1), 'd')
dae.g += spmatrix(
mul(
self.u0,
else:
marker = ''
numeral.append('#' + str(idx + 1) + marker)
# compute frequency, undamped frequency and damping
freq = [0] * neig
ufreq = [0] * neig
damping = [0] * neig
for idx, item in enumerate(mu):
if item.imag == 0:
freq[idx] = 0
ufreq[idx] = 0
damping[idx] = 0
else:
freq[idx] = abs(item) / 2 / pi
ufreq[idx] = abs(item.imag / 2 / pi)
damping[idx] = -div(item.real, abs(item)) * 100
# obtain most associated variables
var_assoc = []
for prow in range(neig):
temp_row = partfact[prow, :]
name_idx = list(temp_row).index(max(temp_row))
var_assoc.append(system.VarName.unamex[name_idx])
pf = []
for prow in range(neig):
temp_row = []
for pcol in range(neig):
temp_row.append(round(partfact[prow, pcol], 5))
pf.append(temp_row)
else:
marker = ''
numeral.append('#' + str(idx + 1) + marker)
# compute frequency, undamped frequency and damping
freq = [0] * neig
ufreq = [0] * neig
damping = [0] * neig
for idx, item in enumerate(mu):
if item.imag == 0:
freq[idx] = 0
ufreq[idx] = 0
damping[idx] = 0
else:
freq[idx] = abs(item) / 2 / pi
ufreq[idx] = abs(item.imag / 2 / pi)
damping[idx] = -div(item.real, abs(item)) * 100
# obtain most associated variables
var_assoc = []
for prow in range(neig):
temp_row = partfact[prow, :]
name_idx = list(temp_row).index(max(temp_row))
var_assoc.append(system.varname.unamex[name_idx])
pf = []
for prow in range(neig):
temp_row = []
for pcol in range(neig):
temp_row.append(round(partfact[prow, pcol], 5))
pf.append(temp_row)
def gycall(self, dae):
dae.add_jac(Gy, mul(self.xmu, 1 - dae.x[self.omega_m]), self.vrd,
self.isq)
dae.add_jac(Gy, -mul(self.xmu, 1 - dae.x[self.omega_m]), self.vrq,
self.isd)
dae.add_jac(Gy, -sin(dae.y[self.a]), self.vsd, self.v)
dae.add_jac(Gy, -mul(dae.y[self.v], cos(dae.y[self.a])), self.vsd,
self.a)
dae.add_jac(Gy, cos(dae.y[self.a]), self.vsq, self.v)
dae.add_jac(Gy, -mul(dae.y[self.v], sin(dae.y[self.a])), self.vsq,
self.a)
dae.add_jac(
Gy,
mul(0.5, self.ngen, pi, self.rho, (self.R)**2, (self.Vwn)**3,
div(1, self.mva_mega), (dae.x[self.vw])**3), self.pw, self.cp)
dae.add_jac(
Gy,
mul(-25.52, (dae.y[self.ilamb])**-2,
exp(mul(-12.5, div(1, dae.y[self.ilamb])))) + mul(
12.5, (dae.y[self.ilamb])**-2,
-1.1 + mul(25.52, div(1, dae.y[self.ilamb])) + mul(
-0.088, dae.x[self.theta_p]),
exp(mul(-12.5, div(1, dae.y[self.ilamb])))), self.cp,
self.ilamb)
dae.add_jac(
Gy,
mul((dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p]))**-2,
(div(1, dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p])) +
mul(-0.035, div(1, 1 + (dae.x[self.theta_p])**3)))**-2),
self.ilamb, self.lamb)
def fcall(self, dae):
toSb = self.Sn / self.system.mva
omega = not0(dae.x[self.omega_m])
dae.f[self.theta_p] = mul(
div(1, self.Tp), -dae.x[self.theta_p] + mul(
self.Kp, self.phi, -1 + dae.x[self.omega_m]))
dae.anti_windup(self.theta_p, 0, pi)
dae.f[self.omega_m] = mul(
0.5, div(1, self.H),
mul(dae.y[self.pw], div(1, dae.x[self.omega_m])) - mul(
self.xmu,
mul(dae.x[self.irq], dae.y[self.isd]) - mul(
dae.x[self.ird], dae.y[self.isq])))
dae.f[self.ird] = mul(
div(1, self.Ts),
-dae.x[self.ird] + mul(self.KV, dae.y[self.v] - dae.y[self.vref]) -
mul(dae.y[self.v], div(1, self.xmu)))
dae.anti_windup(self.ird, self.ird_min, self.irq_max)
k = mul(self.x0, toSb, div(1, dae.y[self.v]), div(1, self.xmu),
div(1, omega))
dae.f[self.irq] = mul(
def gcall(self, dae):
dae.g[self.pw] = -dae.y[self.pw] + mul(
0.5, dae.y[self.cp], self.ngen, pi, self.rho, (self.R)**2,
(self.Vwn)**3, div(1, self.mva_mega), (dae.x[self.vw])**3)
dae.g[self.lamb] = -dae.y[self.lamb] + mul(
4, self.R, self.fn, self.ngb, dae.x[self.omega_m], pi,
div(1, self.Vwn), div(1, self.npole), div(1, dae.x[self.vw]))
dae.g[self.cp] = -dae.y[self.cp] + mul(
-1.1 + mul(25.52, dae.y[self.ilamb]) + mul(-0.08800000000000001,
dae.x[self.theta_p]),
exp(mul(-12.5, dae.y[self.ilamb])))
dae.g[self.ilamb] = div(
1, dae.y[self.lamb] + mul(
0.08, dae.x[self.theta_p])) - dae.y[self.ilamb] + mul(
-0.035, div(1, 1 + (dae.x[self.theta_p])**3))
def wb(self):
"""System base radial frequency"""
return 2 * pi * self.config.freq
Gx,
mul(1.5, dae.y[self.cp],
self.ngen, pi, self.rho, (self.R)**2, (self.Vwn)**3,
div(1, self.mva_mega), (dae.x[self.vw])**2), self.pw, self.vw)
dae.add_jac(
Gx,
mul(-0.088, exp(
mul(-12.5, div(1, dae.y[self.ilamb])))), self.cp, self.theta_p)
dae.add_jac(
Gx,
mul(-4, self.R, self.fn, self.ngb, dae.x[self.omega_m], pi,
div(1, self.Vwn), div(1, self.npole), (dae.x[self.vw])**-2),
self.lamb, self.vw)
dae.add_jac(
Gx,
mul(4, self.R, self.fn, self.ngb, pi, div(1, self.Vwn),
div(1, self.npole), div(1, dae.x[self.vw])), self.lamb,
self.omega_m)
dae.add_jac(
Gx,
mul((div(1, dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p])) +
mul(-0.035, div(1, 1 + (dae.x[self.theta_p]) ** 3))) ** -2,
mul(0.08, (dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p]))
** -2) + mul(-0.105, (dae.x[self.theta_p]) ** 2,
(1 + (dae.x[self.theta_p]) ** 3) ** -2)),
self.ilamb,
self.theta_p)
dae.add_jac(Gx, -mul(self.u0, dae.y[self.vrq]), self.a, self.irq)
dae.add_jac(Gx, -mul(self.u0, dae.y[self.vrd]), self.a, self.ird)
dae.add_jac(Gx, mul(self.u0, dae.y[self.v], self.xmu, div(1, self.x0)),
self.v, self.ird)
def phi(self):
deg1 = pi / 180
dae = self.system.dae
above = agtb(dae.x[self.omega_m], 1)
phi_degree_step = mfloor((dae.x[self.omega_m] - 1) / deg1) * deg1
return mul(phi_degree_step, above)