|
| 1 | +""" |
| 2 | + isfixpoint(nw::Network, s0::NWState; tol=1e-10) |
| 3 | +
|
| 4 | +Check if the state `s0` is a fixpoint of the network `nw` by |
| 5 | +calculating the the RHS and check that every entry is within the |
| 6 | +given tolerance `tol`. |
| 7 | +""" |
| 8 | +function isfixpoint(nw::Network, s0::NWState; tol=1e-10) |
| 9 | + # Check if the state is a fixpoint of the network |
| 10 | + u0 = uflat(s0) |
| 11 | + p0 = pflat(s0) |
| 12 | + du = zeros(eltype(u0), length(u0)) |
| 13 | + nw(du, u0, p0, s0.t) |
| 14 | + |
| 15 | + # Check if the change is within the tolerance |
| 16 | + return all(abs.(du) .< tol) |
| 17 | +end |
| 18 | + |
| 19 | +""" |
| 20 | + is_linear_stable(nw::Network, s0::NWState; kwargs...) |
| 21 | +
|
| 22 | +Check if the fixpoint `s0` of the network `nw` is linearly stable by computing |
| 23 | +the eigenvalues of the Jacobian matrix (or reduced Jacobian for constrained systems). |
| 24 | +
|
| 25 | +A fixpoint is linearly stable if all eigenvalues of the Jacobian have negative |
| 26 | +real parts. For systems with algebraic constraints (non-identity mass matrix), |
| 27 | +the reduced Jacobian is used following the approach in [1]. |
| 28 | +See [`jacobian_eigenvals`](@ref) for more details. |
| 29 | +
|
| 30 | +# Arguments |
| 31 | +- `nw::Network`: The network dynamics object |
| 32 | +- `s0::NWState`: The state to check for linear stability (must be a fixpoint) |
| 33 | +- `kwargs...`: Additional keyword arguments passed to `jacobian_eigenvals` |
| 34 | +
|
| 35 | +# Returns |
| 36 | +- `Bool`: `true` if the fixpoint is linearly stable, `false` otherwise |
| 37 | +
|
| 38 | +# References |
| 39 | +[1] "Power System Modelling and Scripting", F. Milano, Chapter 7.2. |
| 40 | +""" |
| 41 | +function is_linear_stable(nw::Network, s0; kwargs...) |
| 42 | + isfixpoint(nw, s0) || error("The state s0 is not a fixpoint of the network nw.") |
| 43 | + λ = jacobian_eigenvals(nw, s0; kwargs...) |
| 44 | + if all(λ -> real(λ) < 0.0, λ) |
| 45 | + return true |
| 46 | + else |
| 47 | + return false |
| 48 | + end |
| 49 | +end |
| 50 | + |
| 51 | +""" |
| 52 | + jacobian_eigenvals(nw::Network, s0::NWState; eigvalf=LinearAlgebra.eigvals) |
| 53 | +
|
| 54 | +Compute the eigenvalues of the Jacobian matrix for linear stability analysis of |
| 55 | +the network dynamics at state `s0`. |
| 56 | +
|
| 57 | +For systems without algebraic constraints (identity mass matrix), this returns |
| 58 | +the eigenvalues of the full Jacobian matrix. For constrained systems (non-identity |
| 59 | +mass matrix), it computes the eigenvalues of the reduced Jacobian following the |
| 60 | +approach for differential-algebraic equations outlined in [1] |
| 61 | +
|
| 62 | +# Arguments |
| 63 | +- `nw::Network`: The network dynamics object |
| 64 | +- `s0::NWState`: The state at which to compute the Jacobian eigenvalues |
| 65 | +- `eigvalf`: Function to compute eigenvalues (default: `LinearAlgebra.eigvals`) |
| 66 | +
|
| 67 | +# Returns |
| 68 | +- `Vector`: Eigenvalues of the Jacobian (or reduced Jacobian for constrained systems) |
| 69 | +
|
| 70 | +# Algorithm |
| 71 | +For unconstrained systems (M = I): |
| 72 | +- Computes eigenvalues of the full Jacobian J |
| 73 | +
|
| 74 | +For constrained systems (M ≠ I, differential-algebraic equations): |
| 75 | +- The system has the form: M * dz/dt = f(z, t) where M is a diagonal mass matrix |
| 76 | +- Variables are partitioned into differential (M_ii = 1) and algebraic (M_ii = 0) components |
| 77 | +- Let z = [x; y] where x are differential and y are algebraic variables |
| 78 | +- The Jacobian J = ∂f/∂z is partitioned as: |
| 79 | + ``` |
| 80 | + J = [f_x f_y] where f_x = ∂f_d/∂x, f_y = ∂f_d/∂y |
| 81 | + [g_x g_y] g_x = ∂g_a/∂x, g_y = ∂g_a/∂y |
| 82 | + ``` |
| 83 | +- For the algebraic constraints 0 = g_a(x, y), we have dy/dt = -g_y^(-1) * g_x * dx/dt |
| 84 | +- Substituting into the differential equations gives the reduced system: |
| 85 | + dx/dt = (f_x - f_y * g_y^(-1) * g_x) * x = A_s * x |
| 86 | +- The eigenvalues of the reduced Jacobian A_s determine stability |
| 87 | +- This approach follows the theory of differential-algebraic equations [1] |
| 88 | +
|
| 89 | +# References |
| 90 | +[1] "Power System Modelling and Scripting", F. Milano, Chapter 7.2. |
| 91 | +""" |
| 92 | +function jacobian_eigenvals(nw::Network, s0; eigvalf=LinearAlgebra.eigvals) |
| 93 | + x0, p, t = uflat(s0), pflat(s0), s0.t # unpack state |
| 94 | + M = nw.mass_matrix |
| 95 | + |
| 96 | + h!(dx, x) = nw(dx, x, p, t) |
| 97 | + j(x) = (dx = similar(x); ForwardDiff.jacobian(h!, dx, x)) # Jacobian |
| 98 | + J = j(x0) # Full Jacobian at the equilibrium point |
| 99 | + |
| 100 | + if M == LinearAlgebra.I # Constraint free system -> use eigenvalues of jacobian |
| 101 | + return eigvalf(J) |
| 102 | + else # Constraints -> Use eigenvalues of reduced jacobian |
| 103 | + M != LinearAlgebra.Diagonal(M) && error("The constraints are not diagonal.") |
| 104 | + c_idx = findall(LinearAlgebra.diag(M) .== 0) |
| 105 | + d_idx = findall(LinearAlgebra.diag(M) .== 1) |
| 106 | + |
| 107 | + f_x = J[d_idx, d_idx] # Differential equations evaluated at the differential variables |
| 108 | + f_y = J[d_idx, c_idx] # Differential equations evaluated at the constrained variables |
| 109 | + |
| 110 | + g_x = J[c_idx, d_idx] # Constrained equations evaluated at the differential variables |
| 111 | + g_y = J[c_idx, c_idx] # Constrained equations evaluated at the constrained variables |
| 112 | + |
| 113 | + D = f_y * LinearAlgebra.pinv(g_y) * g_x # Degradation matrix |
| 114 | + A_s = f_x - D # State matrix / Reduced Jacobian (eq. 7.16 in [1]) |
| 115 | + return eigvalf(A_s) # Eigenvalues of the reduced jacobian |
| 116 | + end |
| 117 | +end |
0 commit comments