Symbolic Analysis with SymPy#
Information |
Details |
|---|---|
Lead Author |
Hantao Cui |
Learning Objectives |
• Understand symbolic computation in power systems |
Prerequisites |
• Python, and NumPy basics |
Estimated Time |
90-240 minutes |
Topics |
• Symbolic computation basics |
Difficulty |
Advanced (4/5) |
Development Status |
Completed |
Last Updated on |
April 22, 2025 |
Intro to Symbolic Computation in Power Systems#
Power system analysis often involves complex mathematical equations that describe the behavior of electrical networks. While numerical methods are commonly used to solve these equations, symbolic computation offers a more powerful approach for understanding the fundamental relationships between system variables and parameters.
Symbolic computation allows us to work with mathematical expressions in their exact form, maintaining variables as symbols rather than specific numerical values. This approach provides several key advantages for power system engineers:
It reveals the underlying structure of equations, making it easier to understand how different parameters affect system behavior.
It enables analytical derivation of sensitivities and relationships that might be obscured in purely numerical approaches.
It allows for automatic differentiation and simplification of complex expressions.
It can generate optimized code for numerical evaluation, potentially improving computational efficiency.
SymPy is a Python library for symbolic mathematics that provides these capabilities. It enables us to create symbolic variables, perform algebraic manipulations, differentiate expressions, solve equations symbolically, and convert symbolic expressions into efficient numerical code.
In power systems, symbolic computation can be particularly valuable for tasks such as:
Deriving the Jacobian matrix for power flow analysis
Understanding the sensitivity of system states to parameter changes
Deriving analytical expressions for contingency analysis
Simplifying complex network equations to gain insights
This module demonstrates how to integrate SymPy with object-oriented power system models to enhance both understanding and computational capabilities.
Learning Objectives
By the end of this module, you will be able to:
Create symbolic representations of power system components and equations
Derive power flow equations and Jacobian matrices symbolically
Understand how symbolic computation can provide insights into system behavior
Convert symbolic expressions to efficient numerical code
Perform sensitivity analysis using symbolic differentiation
Setting Up the Environment#
Before we begin working with symbolic computation, we need to import the necessary libraries and review the power system component classes we’ll be using.
%run oop-and-dc-power-flow.ipynb
import numpy as np
import sympy as sp
from sympy import symbols, Matrix, cos, sin, lambdify
The power system component classes (Bus, Branch, and Generator) that we developed in the previous module provide an object-oriented representation of the physical components in a power system. The %run magic brought these classes into this notebook, so we won’t need to redefine them.
In this module, we’ll enhance these classes by incorporating symbolic computation to gain deeper insights into system behavior and relationships.
Introduction to SymPy#
Before diving into symbolic analysis of power systems, it is essential to understand the basics of SymPy, the symbolic mathematics library in Python that allows us to work with mathematical expressions in their exact, symbolic form rather than with numerical approximations.
Key SymPy Components#
Creating Symbolic Variables#
import sympy as sp
# Create a single symbolic variable
x = sp.Symbol('x')
y = sp.Symbol('y')
# Create multiple symbolic variables at once
a, b, c = sp.symbols('a b c')
# Create indexed variables
theta = sp.symbols('theta_1:4') # Creates theta_1, theta_2, theta_3
The Symbol and symbols functions create symbolic variables that can be manipulated algebraically. Unlike numerical variables, these maintain their symbolic nature through mathematical operations.
Symbolic Matrices#
# Create a 2x2 symbolic matrix
A = sp.Matrix([[a, b], [c, 0]])
# Matrix operations
A_inverse = A.inv() # Symbolic matrix inverse
det_A = A.det() # Symbolic determinant
print(f"A's determinant: {det_A}")
A's determinant: -b*c
SymPy’s Matrix class allows creation and manipulation of matrices containing symbolic expressions. This is essential for power system analysis, where we work with matrices like the admittance matrix (Y-bus) and Jacobian.
Mathematical Functions#
# Trigonometric functions
expr1 = sp.sin(x) * sp.cos(y)
# Complex numbers
z = a + sp.I * b # I is the imaginary unit
real_part = sp.re(z) # Extract real part
imag_part = sp.im(z) # Extract imaginary part
SymPy provides symbolic versions of mathematical functions that differ from NumPy or Python’s built-in functions. These maintain symbolic expressions rather than computing numerical values. The sp.I represents the imaginary unit (√-1), and functions like re and im extract real and imaginary parts of complex expressions.
Differentiation#
# Differentiate an expression with respect to a variable
f = x**2 * sp.sin(y)
df_dx = sp.diff(f, x) # Partial derivative with respect to x
df_dy = sp.diff(f, y) # Partial derivative with respect to y
print(f"df/dx = {df_dx}")
print(f"df/dy = {df_dy}")
df/dx = 2*x*sin(y)
df/dy = x**2*cos(y)
The diff function computes symbolic derivatives, which is essential for deriving the Jacobian matrix in power flow analysis.
Simplification and Substitution#
# Simplify complex expressions
expr = sp.sin(x)**2 + sp.cos(x)**2
simplified = sp.simplify(expr) # Simplifies to 1
# Substitute values or expressions
expr = a*x + b*y
subs_expr = expr.subs([(x, 2), (y, a)]) # Substitutes x=2 and y=a
print(f"expr = {expr}")
print(f"subs_expr = {subs_expr}")
expr = a*x + b*y
subs_expr = a*b + 2*a
SymPy can simplify expressions and substitute values or other expressions into symbolic formulas. This is useful for applying assumptions (like in DC power flow) or evaluating expressions at specific points.
Converting to Numerical Functions#
# Create a symbolic expression
expr = sp.sin(x) * sp.cos(y) + x**2
# Convert to a numerical function
f_numeric = sp.lambdify([x, y], expr, 'numpy')
# Now can be evaluated with numerical values
result = f_numeric(0.5, 1.0) # Fast numerical evaluation
print(f"expr = {expr}")
print(f"f_numeric = {f_numeric}")
expr = x**2 + sin(x)*cos(y)
f_numeric = <function _lambdifygenerated at 0x7fcc0b0d2020>
print("Source code of `f_numeric` using `inspect`:")
from inspect import getsource
print(getsource(f_numeric))
Source code of `f_numeric` using `inspect`:
def _lambdifygenerated(x, y):
return x**2 + sin(x)*cos(y)
The lambdify function converts symbolic expressions to fast numerical functions using libraries like NumPy. The name lambdify means turning a symbolic expression into a lambda function. This bridges the gap between symbolic analysis and efficient numerical computation, allowing us to derive expressions symbolically but evaluate them quickly.
Differences from NumPy Functions#
A key distinction is that SymPy functions (like sin, cos) maintain expressions symbolically:
import numpy as np
# NumPy computes a numerical value
np_result = np.sin(0.5) # Returns 0.479425...
# SymPy maintains a symbolic expression
x = sp.Symbol('x')
sp_result = sp.sin(x) # Returns sin(x) as a symbolic expression
sp_eval = sp.sin(0.5) # Returns sin(0.5) exactly or as a high-precision approximation
In power system analysis, this symbolic approach helps us:
Derive general formulas that apply to any system
Understand the mathematical structure of power flow equations
Perform analytical sensitivity analysis
Generate optimized code for specific numerical problems
With this foundation in SymPy, we can now apply these techniques to symbolic analysis of power systems, gaining deeper insights while maintaining computational efficiency.
Symbolic Representation of Power System Components#
One of the key advantages of symbolic computation is the ability to work with variables as abstract symbols rather than specific numerical values. This allows us to derive general expressions that hold true for any valid input values.
Let’s start by creating symbolic variables for the key quantities in power system analysis:
# Create symbolic variables for bus voltages and angles
def create_symbolic_bus_variables(bus_count):
# Create voltage magnitude symbols for each bus
V = [symbols(f'V_{i}') for i in range(1, bus_count + 1)]
# Create voltage angle symbols for each bus
theta = [symbols(f'θ_{i}') for i in range(1, bus_count + 1)]
return V, theta
# Create symbolic variables for branch parameters
def create_symbolic_branch_parameters(branch_count):
# Create resistance symbols
R = [symbols(f'R_{i}') for i in range(1, branch_count + 1)]
# Create reactance symbols
X = [symbols(f'X_{i}') for i in range(1, branch_count + 1)]
# Create susceptance symbols
B = [symbols(f'B_{i}') for i in range(1, branch_count + 1)]
return R, X, B
With these symbolic variables, we can extend our power system component classes to work with both numerical and symbolic representations. For example, we can create a symbolic version of the power flow equations:
def symbolic_power_flow_equations(V, theta, Y_bus):
"""
Create symbolic expressions for the active and reactive power flow equations.
Parameters:
-----------
V : list of sympy symbols
Bus voltage magnitudes
theta : list of sympy symbols
Bus voltage angles
Y_bus : Matrix of complex values or symbols
Bus admittance matrix (can be numeric or symbolic)
Returns:
--------
P, Q : lists of symbolic expressions
Active and reactive power injections at each bus
"""
n = len(V)
P = []
Q = []
for i in range(n):
p_i = 0
q_i = 0
for j in range(n):
# Extract magnitude and angle from complex Y_bus
if isinstance(Y_bus[i, j], complex):
G_ij = sp.re(Y_bus[i, j])
B_ij = sp.im(Y_bus[i, j])
else:
# If Y_bus is symbolic
G_ij = sp.symbols(f'G_{i+1}{j+1}')
B_ij = sp.symbols(f'B_{i+1}{j+1}')
# Active power equation term
p_term = V[i] * V[j] * (G_ij * cos(theta[i] - theta[j]) +
B_ij * sin(theta[i] - theta[j]))
# Reactive power equation term
q_term = V[i] * V[j] * (G_ij * sin(theta[i] - theta[j]) -
B_ij * cos(theta[i] - theta[j]))
p_i += p_term
q_i += q_term
P.append(p_i)
Q.append(q_i)
return P, Q
# Example: Creating symbolic variables for a 3-bus system
bus_count = 3
V, theta = create_symbolic_bus_variables(bus_count)
# Display the created symbolic variables
print("Voltage magnitude symbols:")
for v in V:
print(f" {v}")
print("\nVoltage angle symbols:")
for t in theta:
print(f" {t}")
# Create symbolic branch parameters for a system with 2 branches
branch_count = 2
R, X, B = create_symbolic_branch_parameters(branch_count)
# Display the branch parameters
print("\nBranch resistance symbols:")
for r in R:
print(f" {r}")
print("\nBranch reactance symbols:")
for x in X:
print(f" {x}")
Voltage magnitude symbols:
V_1
V_2
V_3
Voltage angle symbols:
θ_1
θ_2
θ_3
Branch resistance symbols:
R_1
R_2
Branch reactance symbols:
X_1
X_2
This function creates symbolic expressions for the active and reactive power injections at each bus, based on the standard AC power flow equations. These expressions maintain the mathematical relationships in symbolic form, allowing us to analyze them analytically.
Symbolic vs. Numeric Representation
When working with power systems, it’s valuable to understand the distinction between symbolic and numerical representations:
Symbolic representation: Maintains variables as abstract symbols, allowing for algebraic manipulation and exact mathematical operations.
Numeric representation: Substitutes specific values for variables, allowing for numerical calculations and specific system solutions.
SymPy allows us to seamlessly transition between these two representations, giving us the best of both worlds: the clarity and insight of symbolic analysis and the computational efficiency of numerical methods.
Admittance Matrix Formation with SymPy#
The admittance matrix (Y-bus) is a fundamental component of power system analysis, representing the network connectivity and impedance characteristics. Let’s see how we can create a symbolic representation of the Y-bus matrix:
# Example: Creating symbolic Y-bus matrix using the PowerSystem class
import numpy as np
import sympy as sp
from sympy import symbols, Matrix, I, pprint
def create_symbolic_y_bus_vectorized(system):
"""
Create a symbolic admittance matrix using vectorized operations.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
Y_bus : Matrix of symbolic expressions
Symbolic bus admittance matrix
"""
n_bus = system.bus.get_count()
Y_bus = sp.zeros(n_bus, n_bus)
# Create a mapping from bus numbers to indices
bus_idx = {int(bus_i): i for i, bus_i in enumerate(system.bus.bus_i)}
# Process all branches at once using vectorized operations
for idx in range(system.branch.get_count()):
# Get from and to bus indices
from_idx = bus_idx[int(system.branch.fbus[idx])]
to_idx = bus_idx[int(system.branch.tbus[idx])]
# Create symbolic impedance parameters
r = sp.symbols(f'r_{from_idx+1}_{to_idx+1}')
x = sp.symbols(f'x_{from_idx+1}_{to_idx+1}')
b = sp.symbols(f'b_{from_idx+1}_{to_idx+1}')
# Calculate admittance
y = 1 / (r + I * x)
# Add to Y-bus matrix
Y_bus[from_idx, to_idx] -= y
Y_bus[to_idx, from_idx] -= y # Maintain symmetry
# Add shunt admittance
Y_bus[from_idx, from_idx] += y + I * b / 2
Y_bus[to_idx, to_idx] += y + I * b / 2
return Y_bus
# Load the system data using the existing parser
system = load_data('case3.xlsx')
# Create symbolic Y-bus matrix
Y_bus_symbolic = create_symbolic_y_bus_vectorized(system)
# Display the symbolic Y-bus matrix
print("Symbolic Y-bus matrix using the PowerSystem class:")
pprint(Y_bus_symbolic)
Symbolic Y-bus matrix using the PowerSystem class:
⎡ⅈ⋅b₁ ₂ ⅈ⋅b₃ ₁ 1 1 -1 ↪
⎢────── + ────── + ───────────── + ───────────── ─────────── ↪
⎢ 2 2 r₃ ₁ + ⅈ⋅x₃ ₁ r₁ ₂ + ⅈ⋅x₁ ₂ r₁ ₂ + ⅈ⋅x₁ ↪
⎢ ↪
⎢ -1 ⅈ⋅b₁ ₂ ⅈ⋅b₂ ₃ 1 ↪
⎢ ───────────── ────── + ────── + ────────── ↪
⎢ r₁ ₂ + ⅈ⋅x₁ ₂ 2 2 r₂ ₃ + ⅈ⋅x ↪
⎢ ↪
⎢ -1 -1 ↪
⎢ ───────────── ─────────── ↪
⎣ r₃ ₁ + ⅈ⋅x₃ ₁ r₂ ₃ + ⅈ⋅x₂ ↪
↪ -1 ⎤
↪ ── ───────────── ⎥
↪ ₂ r₃ ₁ + ⅈ⋅x₃ ₁ ⎥
↪ ⎥
↪ 1 -1 ⎥
↪ ─── + ───────────── ───────────── ⎥
↪ ₂ ₃ r₁ ₂ + ⅈ⋅x₁ ₂ r₂ ₃ + ⅈ⋅x₂ ₃ ⎥
↪ ⎥
↪ ⅈ⋅b₂ ₃ ⅈ⋅b₃ ₁ 1 1 ⎥
↪ ── ────── + ────── + ───────────── + ─────────────⎥
↪ ₃ 2 2 r₃ ₁ + ⅈ⋅x₃ ₁ r₂ ₃ + ⅈ⋅x₂ ₃⎦
This function creates a symbolic representation of the Y-bus matrix, using symbols to represent impedances and admittances. The resulting matrix maintains the mathematical structure while allowing for symbolic manipulation.
We can also create a more specific symbolic Y-bus by incorporating the actual branch parameters:
# Can also create a numerical version with specific parameter values
def create_numerical_y_bus(system):
"""
Create a numerical admittance matrix using the system data.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
Y_bus : numpy array of complex values
Numerical bus admittance matrix
"""
n_bus = system.bus.get_count()
Y_bus = np.zeros((n_bus, n_bus), dtype=complex)
# Create a mapping from bus numbers to indices
bus_idx = {int(bus_i): i for i, bus_i in enumerate(system.bus.bus_i)}
# Process all branches at once
for idx in range(system.branch.get_count()):
# Get from and to bus indices
from_idx = bus_idx[int(system.branch.fbus[idx])]
to_idx = bus_idx[int(system.branch.tbus[idx])]
# Calculate admittance
z = complex(system.branch.r[idx], system.branch.x[idx])
y = 1 / z
# Add to Y-bus matrix
Y_bus[from_idx, to_idx] -= y
Y_bus[to_idx, from_idx] -= y # Maintain symmetry
# Add shunt admittance
Y_bus[from_idx, from_idx] += y + 1j * system.branch.b[idx] / 2
Y_bus[to_idx, to_idx] += y + 1j * system.branch.b[idx] / 2
return Y_bus
# Create numerical Y-bus matrix
Y_bus_numerical = create_numerical_y_bus(system)
# Display the numerical Y-bus matrix
print("\nNumerical Y-bus matrix using the system data:")
print(np.round(Y_bus_numerical, 4)) # Round for cleaner display
Numerical Y-bus matrix using the system data:
[[ 1.2113-13.0946j -0.2212 +3.3186j -0.9901 +9.901j ]
[-0.2212 +3.3186j 1.2113-13.0946j -0.9901 +9.901j ]
[-0.9901 +9.901j -0.9901 +9.901j 1.9802-19.702j ]]
# Now combine symbolic and numerical approaches
def create_symbolic_from_numerical(Y_bus_numerical):
"""
Create a symbolic Y-bus with numerical values substituted.
Parameters:
-----------
Y_bus_numerical : numpy array
Numerical Y-bus matrix
Returns:
--------
Y_bus_symbolic : Matrix of symbolic expressions
Symbolic Y-bus matrix with numerical values
"""
n = Y_bus_numerical.shape[0]
Y_bus_symbolic = sp.zeros(n, n)
for i in range(n):
for j in range(n):
if Y_bus_numerical[i, j] != 0:
# Create symbolic variable for non-zero elements
var_name = f'Y_{i+1}{j+1}'
Y_ij = sp.symbols(var_name)
# Substitute with numerical value
Y_bus_symbolic[i, j] = Y_ij
# Create substitution dictionary
subs_dict = {}
for i in range(n):
for j in range(n):
if Y_bus_numerical[i, j] != 0:
var_name = f'Y_{i+1}{j+1}'
subs_dict[sp.symbols(var_name)] = complex(np.round(Y_bus_numerical[i, j].real, 4),
np.round(Y_bus_numerical[i, j].imag, 4))
return Y_bus_symbolic, subs_dict
# Create symbolic Y-bus with numerical values
Y_bus_sym_with_num, subs_dict = create_symbolic_from_numerical(Y_bus_numerical)
print("\nSymbolic Y-bus matrix with variables:")
pprint(Y_bus_sym_with_num)
print("\nSubstitution dictionary for numerical values:")
for key, value in subs_dict.items():
print(f"{key} = {value}")
Symbolic Y-bus matrix with variables:
⎡Y₁₁ Y₁₂ Y₁₃⎤
⎢ ⎥
⎢Y₂₁ Y₂₂ Y₂₃⎥
⎢ ⎥
⎣Y₃₁ Y₃₂ Y₃₃⎦
Substitution dictionary for numerical values:
Y_11 = (1.2113-13.0946j)
Y_12 = (-0.2212+3.3186j)
Y_13 = (-0.9901+9.901j)
Y_21 = (-0.2212+3.3186j)
Y_22 = (1.2113-13.0946j)
Y_23 = (-0.9901+9.901j)
Y_31 = (-0.9901+9.901j)
Y_32 = (-0.9901+9.901j)
Y_33 = (1.9802-19.702j)
The symbolic Y-bus matrix provides insights into the network structure and can be manipulated symbolically to understand the relationships between network parameters and system behavior.
# Example: Calculate symbolic power flow using the Y-bus
def symbolic_power_injections(system, Y_bus_symbolic):
"""
Calculate symbolic power injections using the Y-bus matrix.
Parameters:
-----------
system : PowerSystem object
The power system data
Y_bus_symbolic : Matrix of symbolic expressions
Symbolic bus admittance matrix
Returns:
--------
P, Q : lists of symbolic expressions
Active and reactive power injections at each bus
"""
n_bus = system.bus.get_count()
# Create symbolic variables for voltages and angles
V = [sp.symbols(f'V_{i+1}') for i in range(n_bus)]
theta = [sp.symbols(f'theta_{i+1}') for i in range(n_bus)]
# Calculate power injections
P = []
Q = []
for i in range(n_bus):
p_i = 0
q_i = 0
for j in range(n_bus):
# Extract symbolic Y-bus elements
y_ij = Y_bus_symbolic[i, j]
if y_ij != 0:
# Express in terms of magnitude and angle
if isinstance(y_ij, complex):
G_ij = float(y_ij.real)
B_ij = float(y_ij.imag)
else:
G_ij = sp.re(y_ij)
B_ij = sp.im(y_ij)
# Calculate power terms
p_term = V[i] * V[j] * (G_ij * sp.cos(theta[i] - theta[j]) +
B_ij * sp.sin(theta[i] - theta[j]))
q_term = V[i] * V[j] * (G_ij * sp.sin(theta[i] - theta[j]) -
B_ij * sp.cos(theta[i] - theta[j]))
p_i += p_term
q_i += q_term
P.append(p_i)
Q.append(q_i)
return P, Q
# Calculate symbolic power injections
P_symbolic, Q_symbolic = symbolic_power_injections(system, Y_bus_symbolic)
# Display the first elements of P and Q to see their structure
print("\nSymbolic active power injection at bus 1:")
pprint(P_symbolic[0])
print("\nSymbolic reactive power injection at bus 1:")
pprint(Q_symbolic[0])
Symbolic active power injection at bus 1:
2 ⎛ im(b₁ ₂) im(b₃ ₁) re(r₃ ₁) - im(x₃ ₁) ↪
V₁ ⋅⎜- ──────── - ──────── + ─────────────────────────────────────────────── + ↪
⎜ 2 2 2 2 ↪
⎝ (re(r₃ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) ↪
↪ re(r₁ ₂) - im(x₁ ₂) ⎞ ⎛ (re(r₁ ₂) ↪
↪ ───────────────────────────────────────────────⎟ + V₁⋅V₂⋅⎜- ─────────────── ↪
↪ 2 2⎟ ⎜ ↪
↪ (re(r₁ ₂) - im(x₁ ₂)) + (re(x₁ ₂) + im(r₁ ₂)) ⎠ ⎝ (re(r₁ ₂) - im( ↪
↪ - im(x₁ ₂))⋅cos(θ₁ - θ₂) (-re(x₁ ₂) - im(r₁ ₂))⋅sin(θ₁ - θ₂) ↪
↪ ──────────────────────────────── - ───────────────────────────────────────── ↪
↪ 2 2 2 ↪
↪ x₁ ₂)) + (re(x₁ ₂) + im(r₁ ₂)) (re(r₁ ₂) - im(x₁ ₂)) + (re(x₁ ₂) + im(r ↪
↪ ⎞ ⎛ (re(r₃ ₁) - im(x₃ ₁))⋅cos(θ₁ - θ₃) ( ↪
↪ ──────⎟ + V₁⋅V₃⋅⎜- ─────────────────────────────────────────────── - ─────── ↪
↪ 2⎟ ⎜ 2 2 ↪
↪ ₁ ₂)) ⎠ ⎝ (re(r₃ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) (re(r₃ ↪
↪ -re(x₃ ₁) - im(r₃ ₁))⋅sin(θ₁ - θ₃) ⎞
↪ ────────────────────────────────────────⎟
↪ 2 2⎟
↪ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) ⎠
Symbolic reactive power injection at bus 1:
2 ⎛ re(b₁ ₂) re(b₃ ₁) -re(x₃ ₁) - im(r₃ ₁) ↪
V₁ ⋅⎜- ──────── - ──────── - ─────────────────────────────────────────────── - ↪
⎜ 2 2 2 2 ↪
⎝ (re(r₃ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) ↪
↪ -re(x₁ ₂) - im(r₁ ₂) ⎞ ⎛ (re(r₁ ₂) ↪
↪ ───────────────────────────────────────────────⎟ + V₁⋅V₂⋅⎜- ─────────────── ↪
↪ 2 2⎟ ⎜ ↪
↪ (re(r₁ ₂) - im(x₁ ₂)) + (re(x₁ ₂) + im(r₁ ₂)) ⎠ ⎝ (re(r₁ ₂) - im( ↪
↪ - im(x₁ ₂))⋅sin(θ₁ - θ₂) (-re(x₁ ₂) - im(r₁ ₂))⋅cos(θ₁ - θ₂) ↪
↪ ──────────────────────────────── + ───────────────────────────────────────── ↪
↪ 2 2 2 ↪
↪ x₁ ₂)) + (re(x₁ ₂) + im(r₁ ₂)) (re(r₁ ₂) - im(x₁ ₂)) + (re(x₁ ₂) + im(r ↪
↪ ⎞ ⎛ (re(r₃ ₁) - im(x₃ ₁))⋅sin(θ₁ - θ₃) ( ↪
↪ ──────⎟ + V₁⋅V₃⋅⎜- ─────────────────────────────────────────────── + ─────── ↪
↪ 2⎟ ⎜ 2 2 ↪
↪ ₁ ₂)) ⎠ ⎝ (re(r₃ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) (re(r₃ ↪
↪ -re(x₃ ₁) - im(r₃ ₁))⋅cos(θ₁ - θ₃) ⎞
↪ ────────────────────────────────────────⎟
↪ 2 2⎟
↪ ₁) - im(x₃ ₁)) + (re(x₃ ₁) + im(r₃ ₁)) ⎠
Symbolic Network Insights
Working with a symbolic Y-bus matrix can provide insights that might be difficult to see in a purely numerical representation:
It reveals the structural zeros in the matrix, showing which buses are not directly connected.
It shows how the network connectivity affects the mathematical formulation of the power flow problem.
It allows for analytical derivation of sensitivity factors, such as how changes in line impedances affect power flows.
Power Flow Equations in Symbolic Form#
Power flow analysis is a fundamental technique in power systems engineering, used to determine the state of a power system under steady-state conditions. The power flow equations relate the power injections at each bus to the bus voltages and network parameters.
With SymPy, we can express these equations symbolically and gain insights into their structure and relationships:
def symbolic_ac_power_flow_equations(system):
"""
Create symbolic AC power flow equations for a given power system.
Parameters:
-----------
system : PowerSystem object
The power system containing bus and branch data
Returns:
--------
P_eq, Q_eq : lists of symbolic expressions
Active and reactive power balance equations
V, theta : lists of symbolic variables
Bus voltage magnitudes and angles
"""
n_bus = system.bus.get_count()
# Create symbolic variables for bus voltages and angles
V = [sp.symbols(f'V_{i+1}') for i in range(n_bus)]
theta = [sp.symbols(f'theta_{i+1}') for i in range(n_bus)]
# Create symbolic Y-bus using the vectorized method
Y_bus = create_symbolic_y_bus_vectorized(system)
# Calculate power injections
P = []
Q = []
for i in range(n_bus):
p_i = 0
q_i = 0
for j in range(n_bus):
# Extract Y-bus components
if Y_bus[i, j] != 0:
if isinstance(Y_bus[i, j], complex):
G_ij = float(Y_bus[i, j].real)
B_ij = float(Y_bus[i, j].imag)
else:
G_ij = sp.re(Y_bus[i, j])
B_ij = sp.im(Y_bus[i, j])
# Calculate power terms
angle_diff = theta[i] - theta[j]
p_term = V[i] * V[j] * (G_ij * sp.cos(angle_diff) + B_ij * sp.sin(angle_diff))
q_term = V[i] * V[j] * (G_ij * sp.sin(angle_diff) - B_ij * sp.cos(angle_diff))
p_i += p_term
q_i += q_term
P.append(p_i)
Q.append(q_i)
# Create power balance equations
P_eq = []
Q_eq = []
# Dictionary to map bus numbers to indices
bus_idx = {int(b_i): i for i, b_i in enumerate(system.bus.bus_i)}
for i in range(n_bus):
# Get bus index
bus_num = system.bus.bus_i[i]
# Get generation at this bus
gen_at_bus = [j for j, g_bus in enumerate(system.gen.bus) if int(g_bus) == int(bus_num)]
P_gen = sum(system.gen.Pg[j] for j in gen_at_bus) if gen_at_bus else 0
Q_gen = sum(system.gen.Qg[j] for j in gen_at_bus) if gen_at_bus else 0
# Get load at this bus
P_load = system.bus.Pd[i]
Q_load = system.bus.Qd[i]
# Create power balance equations
P_eq.append(P_gen - P_load - P[i])
Q_eq.append(Q_gen - Q_load - Q[i])
return P_eq, Q_eq, V, theta
This function creates symbolic expressions for the power balance equations at each bus. These equations represent the conservation of power: the power injected at a bus (from generation) minus the power consumed (by loads) must equal the power flowing out through the connected branches.
We can also derive the Jacobian matrix symbolically, which is essential for the Newton-Raphson method used to solve the power flow equations:
def symbolic_power_flow_jacobian(system, P_eq, Q_eq, V, theta):
"""
Create symbolic Jacobian matrix for power flow equations.
Parameters:
-----------
system : PowerSystem object
The power system data
P_eq, Q_eq : lists of symbolic expressions
Active and reactive power balance equations
V, theta : lists of symbolic variables
Bus voltage magnitudes and angles
Returns:
--------
J : Matrix of symbolic expressions
Symbolic Jacobian matrix
"""
n_bus = system.bus.get_count()
# Identify bus types and indices
slack_idx = next(i for i, type_val in enumerate(system.bus.bus_type) if type_val == 3)
pv_idx = [i for i, type_val in enumerate(system.bus.bus_type) if type_val == 2]
pq_idx = [i for i, type_val in enumerate(system.bus.bus_type) if type_val == 1]
# Extract unknown variables based on bus types
unknown_theta = [theta[i] for i in range(n_bus) if i != slack_idx]
unknown_V = [V[i] for i in pq_idx]
# Extract relevant equations
active_eqs = [P_eq[i] for i in range(n_bus) if i != slack_idx]
reactive_eqs = [Q_eq[i] for i in pq_idx]
# Combine equations and variables
all_eqs = active_eqs + reactive_eqs
all_vars = unknown_theta + unknown_V
# Create Jacobian matrix
J = sp.zeros(len(all_eqs), len(all_vars))
# Fill Jacobian matrix
for i, eq in enumerate(all_eqs):
for j, var in enumerate(all_vars):
J[i, j] = sp.diff(eq, var)
return J
The symbolic Jacobian matrix reveals the sensitivity of the power balance equations to changes in the unknown variables (voltage magnitudes and angles). This can provide insights into the coupling between active and reactive power and the behavior of the power flow solution.
# Example: Creating symbolic power flow equations for a small power system
# Load the system data
system = load_data('case3.xlsx')
# Create symbolic power flow equations
P_eq, Q_eq, V, theta = symbolic_ac_power_flow_equations(system)
# Display the active power equations
print("Symbolic active power equations:")
for i, eq in enumerate(P_eq):
print(f"P_eq[{i+1}] = {eq}")
# Display the reactive power equations
print("\nSymbolic reactive power equations:")
for i, eq in enumerate(Q_eq):
print(f"Q_eq[{i+1}] = {eq}")
# Create symbolic Jacobian matrix
J = symbolic_power_flow_jacobian(system, P_eq, Q_eq, V, theta)
# Display the Jacobian matrix shape
print(f"\nJacobian matrix shape: {J.shape}")
# Display a portion of the Jacobian matrix (first 2x2 block)
print("\nA portion of the symbolic Jacobian matrix:")
for i in range(min(2, J.shape[0])):
for j in range(min(2, J.shape[1])):
print(f"J[{i},{j}] = {J[i,j]}")
# Simplify one element of the Jacobian to see its structure
if J.shape[0] > 0 and J.shape[1] > 0:
simple_elem = sp.simplify(J[0, 0])
print("\nSimplified first element of the Jacobian:")
print(simple_elem)
Symbolic active power equations:
P_eq[1] = -V_1**2*(-im(b_1_2)/2 - im(b_3_1)/2 + (re(r_3_1) - im(x_3_1))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (re(r_1_2) - im(x_1_2))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_1*V_2*(-(re(r_1_2) - im(x_1_2))*cos(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*sin(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_1*V_3*(-(re(r_3_1) - im(x_3_1))*cos(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*sin(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))
P_eq[2] = -V_1*V_2*(-(re(r_1_2) - im(x_1_2))*cos(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*sin(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_2**2*(-im(b_1_2)/2 - im(b_2_3)/2 + (re(r_2_3) - im(x_2_3))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (re(r_1_2) - im(x_1_2))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_2*V_3*(-(re(r_2_3) - im(x_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + 50
P_eq[3] = -V_1*V_3*(-(re(r_3_1) - im(x_3_1))*cos(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*sin(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) - V_2*V_3*(-(re(r_2_3) - im(x_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) - V_3**2*(-im(b_2_3)/2 - im(b_3_1)/2 + (re(r_3_1) - im(x_3_1))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (re(r_2_3) - im(x_2_3))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) - 120
Symbolic reactive power equations:
Q_eq[1] = -V_1**2*(-re(b_1_2)/2 - re(b_3_1)/2 - (-re(x_3_1) - im(r_3_1))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_1_2) - im(r_1_2))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_1*V_2*(-(re(r_1_2) - im(x_1_2))*sin(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*cos(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_1*V_3*(-(re(r_3_1) - im(x_3_1))*sin(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*cos(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))
Q_eq[2] = -V_1*V_2*((re(r_1_2) - im(x_1_2))*sin(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*cos(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_2**2*(-re(b_1_2)/2 - re(b_2_3)/2 - (-re(x_2_3) - im(r_2_3))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_1_2) - im(r_1_2))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_2*V_3*(-(re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
Q_eq[3] = -V_1*V_3*((re(r_3_1) - im(x_3_1))*sin(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*cos(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) - V_2*V_3*((re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) - V_3**2*(-re(b_2_3)/2 - re(b_3_1)/2 - (-re(x_3_1) - im(r_3_1))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_2_3) - im(r_2_3))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) - 50
Jacobian matrix shape: (3, 3)
A portion of the symbolic Jacobian matrix:
J[0,0] = -V_1*V_2*(-(re(r_1_2) - im(x_1_2))*sin(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*cos(theta_1 - theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) - V_2*V_3*((re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
J[0,1] = -V_2*V_3*(-(re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
J[1,0] = -V_2*V_3*((re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
J[1,1] = -V_1*V_3*(-(re(r_3_1) - im(x_3_1))*sin(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*cos(theta_1 - theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) - V_2*V_3*(-(re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*cos(theta_2 - theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
Simplified first element of the Jacobian:
V_2*(V_1*((re(r_1_2) - im(x_1_2))*sin(theta_1 - theta_2) - (re(x_1_2) + im(r_1_2))*cos(theta_1 - theta_2))*((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - V_3*((re(r_2_3) - im(x_2_3))*sin(theta_2 - theta_3) + (re(x_2_3) + im(r_2_3))*cos(theta_2 - theta_3))*((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))/(((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)*((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))
DC Power Flow with Symbolic Variables#
DC power flow is a simplified version of the full AC power flow, making several assumptions to linearize the equations. With symbolic computation, we can derive the DC power flow equations directly from the AC equations, providing insights into the approximations involved:
def symbolic_dc_power_flow_derivation(system):
"""
Derive DC power flow equations symbolically from AC equations.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
P_dc : list of symbolic expressions
DC power flow equations
"""
n_bus = system.bus.get_count()
# Create symbolic variables
V = [sp.symbols(f'V_{i+1}') for i in range(n_bus)]
theta = [sp.symbols(f'theta_{i+1}') for i in range(n_bus)]
# Create symbolic Y-bus
Y_bus = create_symbolic_y_bus_vectorized(system)
# Calculate AC power injections
P_ac = []
for i in range(n_bus):
p_i = 0
for j in range(n_bus):
if Y_bus[i, j] != 0:
if isinstance(Y_bus[i, j], complex):
G_ij = float(Y_bus[i, j].real)
B_ij = float(Y_bus[i, j].imag)
else:
G_ij = sp.re(Y_bus[i, j])
B_ij = sp.im(Y_bus[i, j])
angle_diff = theta[i] - theta[j]
p_term = V[i] * V[j] * (G_ij * sp.cos(angle_diff) + B_ij * sp.sin(angle_diff))
p_i += p_term
P_ac.append(p_i)
# Apply DC power flow assumptions
# 1. V_i = V_j = 1.0 (flat voltage profile)
# 2. sin(theta_i - theta_j) ≈ (theta_i - theta_j)
# 3. cos(theta_i - theta_j) ≈ 1.0
# 4. G_ij << B_ij (ignore conductance)
P_dc = []
for i in range(n_bus):
# Start with AC active power equation
p_i = P_ac[i]
# Apply assumption 1: V_i = V_j = 1.0
p_i = p_i.subs([(v, 1.0) for v in V])
# Apply assumption 3: cos(theta_i - theta_j) ≈ 1.0
p_i = p_i.subs([(sp.cos(theta[i] - theta[j]), 1.0) for j in range(n_bus)])
# Apply assumption 2: sin(theta_i - theta_j) ≈ (theta_i - theta_j)
p_i = p_i.subs([(sp.sin(theta[i] - theta[j]), theta[i] - theta[j]) for j in range(n_bus)])
# Apply assumption 4: Ignore conductance (G_ij = 0)
if isinstance(p_i, sp.Expr):
p_i = p_i.subs([(sp.re(Y_bus[i, j]), 0) for j in range(n_bus) if Y_bus[i, j] != 0])
P_dc.append(sp.simplify(p_i))
return P_dc
This function symbolically derives the DC power flow equations by applying the standard DC approximations to the full AC equations. The resulting expressions show how the DC power flow relates to the angle differences and network susceptances, providing a clear understanding of the linearization process.
We can also create the DC power flow matrix (B’) symbolically:
def symbolic_dc_power_flow_matrix(system):
"""
Create symbolic DC power flow matrix (B').
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
B_prime : Matrix of symbolic expressions
Symbolic DC power flow matrix
"""
n_bus = system.bus.get_count()
B_prime = sp.zeros(n_bus, n_bus)
# Dictionary to map bus numbers to indices
bus_idx = {int(b_i): i for i, b_i in enumerate(system.bus.bus_i)}
# Add branch susceptances to B'
for idx in range(system.branch.get_count()):
fbus = int(system.branch.fbus[idx])
tbus = int(system.branch.tbus[idx])
i = bus_idx[fbus]
j = bus_idx[tbus]
# Create symbolic reactance
x_ij = sp.symbols(f'x_{i+1}_{j+1}')
# Calculate susceptance
b_ij = 1 / x_ij
# Add off-diagonal elements
B_prime[i, j] -= b_ij
B_prime[j, i] -= b_ij
# Add diagonal elements
B_prime[i, i] += b_ij
B_prime[j, j] += b_ij
return B_prime
# Example: Creating DC power flow equations for the test system
# Load the system data
system = load_data('case3.xlsx')
# Create symbolic DC power flow equations
P_dc = symbolic_dc_power_flow_derivation(system)
# Display the DC power flow equations
print("Symbolic DC power flow equations:")
for i, eq in enumerate(P_dc):
print(f"P_dc[{i+1}] = {eq}")
# Create symbolic DC power flow matrix (B')
B_prime = symbolic_dc_power_flow_matrix(system)
# Display the DC power flow matrix
print("\nSymbolic DC power flow matrix (B'):")
for i in range(B_prime.shape[0]):
for j in range(B_prime.shape[1]):
if B_prime[i, j] != 0:
print(f"B'[{i+1},{j+1}] = {B_prime[i, j]}")
# Substitute actual reactance values to get a numerical B' matrix
B_prime_numeric = B_prime.copy()
for idx in range(system.branch.get_count()):
fbus = int(system.branch.fbus[idx])
tbus = int(system.branch.tbus[idx])
i = next(i for i, b_i in enumerate(system.bus.bus_i) if int(b_i) == fbus)
j = next(j for j, b_j in enumerate(system.bus.bus_i) if int(b_j) == tbus)
# Substitute reactance value
x_ij = sp.symbols(f'x_{i+1}_{j+1}')
B_prime_numeric = B_prime_numeric.subs(x_ij, float(system.branch.x[idx]))
print("\nB' matrix with numerical values:")
for i in range(B_prime_numeric.shape[0]):
print([float(B_prime_numeric[i, j]) if B_prime_numeric[i, j].is_number else 0
for j in range(B_prime_numeric.shape[1])])
# Compare with direct calculation of DC power flow
# For a 3-bus system, let's create a simple example with angles
theta_vals = [0.0, 0.1, -0.05] # Example angle values
theta_sym = [sp.symbols(f'theta_{i+1}') for i in range(len(theta_vals))]
# Calculate P_dc with these angle values
P_dc_numeric = []
for eq in P_dc:
# Substitute angle values
eq_subs = eq
for i, theta_val in enumerate(theta_vals):
eq_subs = eq_subs.subs(theta_sym[i], theta_val)
# Try to evaluate to a float
try:
P_dc_numeric.append(float(eq_subs))
except:
P_dc_numeric.append(eq_subs)
print("\nDC power flow results with angle values:")
for i, p in enumerate(P_dc_numeric):
print(f"P_dc[{i+1}] = {p}")
Symbolic DC power flow equations:
P_dc[1] = (-0.5*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0*(im(b_1_2) + im(b_3_1)) + 1.0*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((-re(r_3_1) + 1.0*im(x_3_1))*cos(theta_1 - 1.0*theta_3) + (re(x_3_1) + im(r_3_1))*sin(theta_1 - 1.0*theta_3) + re(r_3_1) - im(x_3_1)) + 1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0*((-re(r_1_2) + 1.0*im(x_1_2))*cos(theta_1 - 1.0*theta_2) + (re(x_1_2) + im(r_1_2))*sin(theta_1 - 1.0*theta_2) + re(r_1_2) - im(x_1_2)))/(((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0)
P_dc[2] = (-0.5*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0*(im(b_1_2) + im(b_2_3)) + 1.0*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((-re(r_2_3) + 1.0*im(x_2_3))*cos(theta_2 - 1.0*theta_3) + (re(x_2_3) + im(r_2_3))*sin(theta_2 - 1.0*theta_3) + re(r_2_3) - im(x_2_3)) - 1.0*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0*(re(x_1_2) + im(r_1_2))*sin(theta_1 - 1.0*theta_2))/(((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0)
P_dc[3] = -1.0*sin(theta_2 - 1.0*theta_3)*re(x_2_3)/((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0 - 1.0*sin(theta_2 - 1.0*theta_3)*im(r_2_3)/((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0 - 1.0*sin(theta_1 - 1.0*theta_3)*re(x_3_1)/((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0 - 1.0*sin(theta_1 - 1.0*theta_3)*im(r_3_1)/((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0 - 0.5*im(b_2_3) - 0.5*im(b_3_1)
Symbolic DC power flow matrix (B'):
B'[1,1] = 1/x_3_1 + 1/x_1_2
B'[1,2] = -1/x_1_2
B'[1,3] = -1/x_3_1
B'[2,1] = -1/x_1_2
B'[2,2] = 1/x_2_3 + 1/x_1_2
B'[2,3] = -1/x_2_3
B'[3,1] = -1/x_3_1
B'[3,2] = -1/x_2_3
B'[3,3] = 1/x_3_1 + 1/x_2_3
B' matrix with numerical values:
[13.333333333333334, -3.3333333333333335, -10.0]
[-3.3333333333333335, 13.333333333333334, -10.0]
[-10.0, -10.0, 20.0]
DC power flow results with angle values:
P_dc[1] = (-0.5*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0*(im(b_1_2) + im(b_3_1)) + 1.0*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*(0.00124973960503372*re(r_3_1) + 0.0499791692706783*re(x_3_1) + 0.0499791692706783*im(r_3_1) - 0.00124973960503372*im(x_3_1)) + 1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0*(0.00499583472197429*re(r_1_2) - 0.0998334166468282*re(x_1_2) - 0.0998334166468282*im(r_1_2) - 0.00499583472197429*im(x_1_2)))/(((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0)
P_dc[2] = (-0.5*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0*(im(b_1_2) + im(b_2_3)) + 1.0*((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*(0.0112289220639578*re(r_2_3) + 0.149438132473599*re(x_2_3) + 0.149438132473599*im(r_2_3) - 0.0112289220639578*im(x_2_3)) + 0.0998334166468282*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0*(re(x_1_2) + im(r_1_2)))/(((re(r_1_2) - 1.0*im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)**1.0*((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0)
P_dc[3] = -0.149438132473599*re(x_2_3)/((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0 - 0.149438132473599*im(r_2_3)/((re(r_2_3) - 1.0*im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)**1.0 - 0.0499791692706783*re(x_3_1)/((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0 - 0.0499791692706783*im(r_3_1)/((re(r_3_1) - 1.0*im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)**1.0 - 0.5*im(b_2_3) - 0.5*im(b_3_1)
The symbolic DC power flow matrix highlights the network structure and the role of branch reactances in determining power flows. It can be used to analytically derive relationships such as Power Transfer Distribution Factors (PTDFs).
DC Power Flow Advantages
DC power flow offers several advantages:
The equations are linear, making them easier to solve
The solution process doesn’t require iteration
The simplicity makes it useful for applications like contingency analysis and market simulations
However, it also has limitations:
It neglects reactive power flows
It assumes all voltages are 1.0 per unit
It neglects transmission losses
Understanding these trade-offs is easier with symbolic analysis, which clearly shows the simplifications made.
Symbolic Sensitivity Analysis#
Sensitivity analysis is crucial for understanding how changes in system parameters or operating conditions affect system behavior. Symbolic computation provides a powerful tool for deriving sensitivity factors analytically.
One important set of sensitivity factors is the Power Transfer Distribution Factors (PTDFs), which show how a change in power injection affects the flow on transmission lines:
def symbolic_ptdf_matrix(system):
"""
Calculate Power Transfer Distribution Factors (PTDFs) symbolically.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
ptdf : Matrix of symbolic expressions
Symbolic PTDF matrix
"""
n_bus = system.bus.get_count()
n_branch = system.branch.get_count()
# Find slack bus index
slack_idx = next(i for i, type_val in enumerate(system.bus.bus_type) if type_val == 3)
# Create symbolic B' matrix
B_prime = symbolic_dc_power_flow_matrix(system)
# Create reduced B' matrix by removing slack bus row and column
B_prime_reduced = B_prime.copy()
B_prime_reduced.row_del(slack_idx)
B_prime_reduced.col_del(slack_idx)
# Calculate B' inverse symbolically
# Note: This can be computationally intensive for large systems
try:
B_prime_inv = B_prime_reduced.inv()
except:
print("Symbolic inversion failed. Using a simplified approach.")
# Use a simplified approach for demonstration
B_prime_inv = sp.zeros(n_bus-1, n_bus-1)
for i in range(n_bus-1):
for j in range(n_bus-1):
B_prime_inv[i, j] = sp.symbols(f'X_{i+1}_{j+1}')
# Create symbolic PTDF matrix
ptdf = sp.zeros(n_branch, n_bus)
# Dictionary to map bus numbers to indices
bus_idx = {int(b_i): i for i, b_i in enumerate(system.bus.bus_i)}
# Calculate PTDF for each branch and bus pair
for l in range(n_branch):
fbus = int(system.branch.fbus[l])
tbus = int(system.branch.tbus[l])
i = bus_idx[fbus]
j = bus_idx[tbus]
# Create symbolic reactance
x_l = sp.symbols(f'x_{i+1}_{j+1}')
for k in range(n_bus):
if k == slack_idx:
# PTDF for slack bus is zero
ptdf[l, k] = 0
else:
# Adjust index for reduced matrix
k_adj = k if k < slack_idx else k - 1
if i != slack_idx and j != slack_idx:
# Both buses are not slack
i_adj = i if i < slack_idx else i - 1
j_adj = j if j < slack_idx else j - 1
ptdf[l, k] = (1 / x_l) * (B_prime_inv[i_adj, k_adj] - B_prime_inv[j_adj, k_adj])
elif i == slack_idx:
# From bus is slack
j_adj = j if j < slack_idx else j - 1
ptdf[l, k] = -(1 / x_l) * B_prime_inv[j_adj, k_adj]
elif j == slack_idx:
# To bus is slack
i_adj = i if i < slack_idx else i - 1
ptdf[l, k] = (1 / x_l) * B_prime_inv[i_adj, k_adj]
return ptdf
This function symbolically derives the PTDF matrix, which relates changes in bus power injections to changes in branch flows. The symbolic derivation reveals how the network structure and impedances determine these sensitivity factors.
We can also derive loss factors symbolically, showing how system losses are affected by changes in bus injections:
def symbolic_loss_factors(system):
"""
Derive loss factors symbolically.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
loss_factors : list of symbolic expressions
Symbolic loss factors for each bus
"""
n_bus = system.bus.get_count()
n_branch = system.branch.get_count()
# Create symbolic variables
V = [sp.symbols(f'V_{i+1}') for i in range(n_bus)]
theta = [sp.symbols(f'theta_{i+1}') for i in range(n_bus)]
# Dictionary to map bus numbers to indices
bus_idx = {int(b_i): i for i, b_i in enumerate(system.bus.bus_i)}
# Calculate system losses symbolically
losses = 0
for idx in range(n_branch):
fbus = int(system.branch.fbus[idx])
tbus = int(system.branch.tbus[idx])
i = bus_idx[fbus]
j = bus_idx[tbus]
# Create symbolic resistance and reactance
r_ij = sp.symbols(f'r_{i+1}_{j+1}')
x_ij = sp.symbols(f'x_{i+1}_{j+1}')
# Calculate impedance and admittance
z_ij = r_ij + sp.I * x_ij
y_ij = 1 / z_ij
# Calculate voltage difference
V_i = V[i] * sp.exp(sp.I * theta[i])
V_j = V[j] * sp.exp(sp.I * theta[j])
# Calculate branch current
I_ij = (V_i - V_j) * y_ij
# Calculate branch losses
branch_loss = r_ij * sp.Abs(I_ij)**2
losses += branch_loss
# Derive loss factors (partial derivative of losses with respect to bus injections)
loss_factors = []
for i in range(n_bus):
# Create symbolic bus injection
P_i = sp.symbols(f'P_{i+1}')
# Loss factor is the partial derivative of losses with respect to P_i
# For simplicity, we'll calculate the partial derivative with respect to theta
loss_factor = 0
# Apply chain rule: ∂Loss/∂P_i = Σj (∂Loss/∂θj) * (∂θj/∂P_i)
for j in range(n_bus):
if j != next(i for i, type_val in enumerate(system.bus.bus_type) if type_val == 3):
# Derivative of losses with respect to theta_j
dlosses_dtheta_j = sp.diff(losses, theta[j])
# Sensitivity of theta_j to change in P_i (simplified as symbolic)
dtheta_j_dP_i = sp.symbols(f'S_{j+1}_{i+1}')
loss_factor += dlosses_dtheta_j * dtheta_j_dP_i
loss_factors.append(loss_factor)
return loss_factors
# Example: Calculating PTDF and loss factors symbolically
# Load the system data
system = load_data('case3.xlsx')
# Calculate symbolic PTDF matrix
try:
ptdf = symbolic_ptdf_matrix(system)
# Display a portion of the PTDF matrix
print("Symbolic PTDF matrix (partial):")
for l in range(min(2, ptdf.shape[0])):
for k in range(min(2, ptdf.shape[1])):
print(f"PTDF[{l+1},{k+1}] = {ptdf[l,k]}")
except Exception as e:
print(f"Could not calculate full symbolic PTDF matrix: {e}")
# Alternative: Calculate PTDF for a single branch-bus pair
print("\nSymbolic PTDF for a single branch-bus pair:")
# For a 3-bus system with bus 1 as slack
# PTDF for branch 1-2 and injection at bus 3
x_12 = sp.symbols('x_1_2')
X_22 = sp.symbols('X_2_2') # Element of B' inverse
ptdf_13_2 = -(1 / x_12) * X_22
print(f"PTDF[branch 1-2, bus 3] = {ptdf_13_2}")
Symbolic PTDF matrix (partial):
PTDF[1,1] = 0
PTDF[1,2] = -(x_1_2*x_2_3 + x_1_2*x_3_1)/(x_1_2*(x_1_2 + x_2_3 + x_3_1))
PTDF[2,1] = 0
PTDF[2,2] = (-x_1_2*x_3_1/(x_1_2 + x_2_3 + x_3_1) + (x_1_2*x_2_3 + x_1_2*x_3_1)/(x_1_2 + x_2_3 + x_3_1))/x_2_3
# Calculate symbolic loss factors
loss_factors = symbolic_loss_factors(system)
# Display the symbolic loss factors
print("\nSymbolic loss factors:")
for i, lf in enumerate(loss_factors):
print(f"Loss Factor[{i+1}] = {lf}")
Symbolic loss factors:
Loss Factor[1] = S_2_1*(2*r_1_2*(r_1_2 + I*x_1_2)*(((re(r_1_2) - im(x_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) + ((re(r_1_2) - im(x_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)))*Abs((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))*sign((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))/(V_1*exp(I*theta_1) - V_2*exp(I*theta_2)) + 2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3))) + S_3_1*(2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3)) + 2*r_3_1*(r_3_1 + I*x_3_1)*(((re(r_3_1) - im(x_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) + ((re(r_3_1) - im(x_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)))*Abs((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))*sign((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))/(V_1*exp(I*theta_1) - V_3*exp(I*theta_3)))
Loss Factor[2] = S_2_2*(2*r_1_2*(r_1_2 + I*x_1_2)*(((re(r_1_2) - im(x_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) + ((re(r_1_2) - im(x_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)))*Abs((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))*sign((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))/(V_1*exp(I*theta_1) - V_2*exp(I*theta_2)) + 2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3))) + S_3_2*(2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3)) + 2*r_3_1*(r_3_1 + I*x_3_1)*(((re(r_3_1) - im(x_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) + ((re(r_3_1) - im(x_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)))*Abs((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))*sign((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))/(V_1*exp(I*theta_1) - V_3*exp(I*theta_3)))
Loss Factor[3] = S_2_3*(2*r_1_2*(r_1_2 + I*x_1_2)*(((re(r_1_2) - im(x_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)) + ((re(r_1_2) - im(x_1_2))*(im(V_1*exp(I*theta_1)) - im(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) + (-re(x_1_2) - im(r_1_2))*(re(V_1*exp(I*theta_1)) - re(V_2*exp(I*theta_2)))/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2))*(-(re(r_1_2) - im(x_1_2))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2) - (-re(x_1_2) - im(r_1_2))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_1_2) - im(x_1_2))**2 + (re(x_1_2) + im(r_1_2))**2)))*Abs((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))*sign((V_1*exp(I*theta_1) - V_2*exp(I*theta_2))/(r_1_2 + I*x_1_2))/(V_1*exp(I*theta_1) - V_2*exp(I*theta_2)) + 2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*((re(r_2_3) - im(x_2_3))*Derivative(im(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(re(V_2*exp(I*theta_2)), theta_2)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3))) + S_3_3*(2*r_2_3*(r_2_3 + I*x_2_3)*(((re(r_2_3) - im(x_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)) + ((re(r_2_3) - im(x_2_3))*(im(V_2*exp(I*theta_2)) - im(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) + (-re(x_2_3) - im(r_2_3))*(re(V_2*exp(I*theta_2)) - re(V_3*exp(I*theta_3)))/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2))*(-(re(r_2_3) - im(x_2_3))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2) - (-re(x_2_3) - im(r_2_3))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_2_3) - im(x_2_3))**2 + (re(x_2_3) + im(r_2_3))**2)))*Abs((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))*sign((V_2*exp(I*theta_2) - V_3*exp(I*theta_3))/(r_2_3 + I*x_2_3))/(V_2*exp(I*theta_2) - V_3*exp(I*theta_3)) + 2*r_3_1*(r_3_1 + I*x_3_1)*(((re(r_3_1) - im(x_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)) + ((re(r_3_1) - im(x_3_1))*(im(V_1*exp(I*theta_1)) - im(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) + (-re(x_3_1) - im(r_3_1))*(re(V_1*exp(I*theta_1)) - re(V_3*exp(I*theta_3)))/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2))*(-(re(r_3_1) - im(x_3_1))*Derivative(im(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2) - (-re(x_3_1) - im(r_3_1))*Derivative(re(V_3*exp(I*theta_3)), theta_3)/((re(r_3_1) - im(x_3_1))**2 + (re(x_3_1) + im(r_3_1))**2)))*Abs((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))*sign((V_1*exp(I*theta_1) - V_3*exp(I*theta_3))/(r_3_1 + I*x_3_1))/(V_1*exp(I*theta_1) - V_3*exp(I*theta_3)))
# Example: Simplify loss factors with assumptions
print("\nSimplified loss factors (partial derivatives only):")
n_bus = system.bus.get_count()
V = [sp.symbols(f'V_{i+1}') for i in range(n_bus)]
theta = [sp.symbols(f'theta_{i+1}') for i in range(n_bus)]
# Create a simplified system loss expression
# L = r_12 * [(V_1*V_2/x_12) * (θ_1 - θ_2)]^2
r_12 = sp.symbols('r_1_2')
x_12 = sp.symbols('x_1_2')
simplified_loss = r_12 * ((V[0]*V[1]/x_12) * (theta[0] - theta[1]))**2
# Calculate partial derivatives
dL_dtheta1 = sp.diff(simplified_loss, theta[0])
dL_dtheta2 = sp.diff(simplified_loss, theta[1])
print(f"∂Loss/∂θ_1 = {dL_dtheta1}")
print(f"∂Loss/∂θ_2 = {dL_dtheta2}")
Simplified loss factors (partial derivatives only):
∂Loss/∂θ_1 = V_1**2*V_2**2*r_1_2*(2*theta_1 - 2*theta_2)/x_1_2**2
∂Loss/∂θ_2 = V_1**2*V_2**2*r_1_2*(-2*theta_1 + 2*theta_2)/x_1_2**2
These symbolic sensitivity factors provide insights into system behavior that might be difficult to obtain through purely numerical methods. They reveal the structural relationships between system parameters and sensitivities.
Practical Considerations
Symbolic computation can be computationally intensive for large systems. For practical applications, it’s often useful to:
Use symbolic computation to derive general formulas
Substitute specific numerical values for final calculations
Generate efficient numerical code from symbolic expressions
Apply simplifications and approximations where appropriate
This approach combines the insight of symbolic analysis with the efficiency of numerical methods.
Converting Symbolic Expressions to Efficient Code#
While symbolic computation provides valuable insights, numerical calculations are often more efficient for practical applications. SymPy provides tools to convert symbolic expressions into efficient numerical functions:
def generate_numerical_functions(symbolic_expressions, symbolic_vars):
"""
Generate efficient numerical functions from symbolic expressions.
Parameters:
-----------
symbolic_expressions : list or Matrix of symbolic expressions
The symbolic expressions to convert
symbolic_vars : list of symbolic variables
The variables in the expressions
Returns:
--------
numeric_functions : list of callable functions
Numerical functions that evaluate the expressions
"""
# Convert a single expression to a list for consistent handling
if not isinstance(symbolic_expressions, list) and not isinstance(symbolic_expressions, Matrix):
symbolic_expressions = [symbolic_expressions]
# Create numerical functions using lambdify
numeric_functions = []
for expr in symbolic_expressions:
# lambdify converts a symbolic expression to a numerical function
numeric_fn = lambdify(symbolic_vars, expr, modules='numpy')
numeric_functions.append(numeric_fn)
return numeric_functions
This function uses SymPy’s lambdify to convert symbolic expressions into efficient numerical functions. These functions can be used for fast repeated evaluations, such as in iterative algorithms like Newton-Raphson.
We can also use this approach to generate optimized code for power system calculations:
def generate_power_flow_functions(system):
"""
Generate efficient numerical functions for power flow calculations.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
power_functions : dict
Dictionary of numerical functions for power calculations
jacobian_function : callable
Function to evaluate the Jacobian matrix
"""
# Create symbolic power flow equations
P_eq, Q_eq, V, theta = symbolic_ac_power_flow_equations(system)
# Create symbolic Jacobian
J = symbolic_power_flow_jacobian(system, P_eq, Q_eq, V, theta)
# Generate numerical functions
all_vars = V + theta
# Convert active power equations to numerical functions
P_functions = []
for eq in P_eq:
# lambdify converts a symbolic expression to a numerical function
P_functions.append(lambdify(all_vars, eq, modules='numpy'))
# Convert reactive power equations to numerical functions
Q_functions = []
for eq in Q_eq:
Q_functions.append(lambdify(all_vars, eq, modules='numpy'))
# Convert Jacobian to numerical function
jacobian_function = lambdify(all_vars, J, modules='numpy')
# Return all functions
power_functions = {
'P': P_functions,
'Q': Q_functions
}
return power_functions, jacobian_function
This function generates optimized numerical functions for power flow calculations, including the power balance equations and the Jacobian matrix. These functions can be used for efficient iterative solution of the power flow problem.
Let’s compare the performance of symbolic-generated code with hard-coded implementations:
def performance_comparison(system):
"""
Compare performance of symbolic-generated and hard-coded implementations.
Parameters:
-----------
system : PowerSystem object
The power system data
Returns:
--------
None (prints comparison results)
"""
import time
import numpy as np
n_bus = system.bus.get_count()
try:
# Generate functions from symbolic expressions
power_functions, jacobian_function = generate_power_flow_functions(system)
# Create test data
V_test = np.ones(n_bus)
theta_test = np.zeros(n_bus)
all_vars_test = np.concatenate([V_test, theta_test])
# Hard-coded implementation for comparison
def hard_coded_power_flow(V, theta, Y_bus):
n = len(V)
P = np.zeros(n)
Q = np.zeros(n)
for i in range(n):
for j in range(n):
if Y_bus[i, j] != 0:
# Extract Y_bus components
G_ij = np.real(Y_bus[i, j])
B_ij = np.imag(Y_bus[i, j])
# Angle difference
angle_diff = theta[i] - theta[j]
# Power calculations
P[i] += V[i] * V[j] * (G_ij * np.cos(angle_diff) +
B_ij * np.sin(angle_diff))
Q[i] += V[i] * V[j] * (G_ij * np.sin(angle_diff) -
B_ij * np.cos(angle_diff))
return P, Q
# Create numerical Y-bus
Y_bus_numerical = np.zeros((n_bus, n_bus), dtype=complex)
# Map bus numbers to indices
bus_idx = {int(b_i): i for i, b_i in enumerate(system.bus.bus_i)}
for idx in range(system.branch.get_count()):
fbus = int(system.branch.fbus[idx])
tbus = int(system.branch.tbus[idx])
i = bus_idx[fbus]
j = bus_idx[tbus]
# Calculate impedance and admittance
r = float(system.branch.r[idx])
x = float(system.branch.x[idx])
b = float(system.branch.b[idx])
z = complex(r, x)
y = 1 / z
# Add off-diagonal elements
Y_bus_numerical[i, j] -= y
Y_bus_numerical[j, i] -= y
# Add diagonal elements
Y_bus_numerical[i, i] += y + 0.5j * b
Y_bus_numerical[j, j] += y + 0.5j * b
# Number of iterations for timing
n_iterations = 100
# Timing for symbolic-generated implementation
start_time = time.time()
for _ in range(n_iterations):
for i in range(n_bus):
if i < len(power_functions['P']):
P_i = power_functions['P'][i](*all_vars_test)
if i < len(power_functions['Q']):
Q_i = power_functions['Q'][i](*all_vars_test)
symbolic_time = time.time() - start_time
# Timing for hard-coded implementation
start_time = time.time()
for _ in range(n_iterations):
P, Q = hard_coded_power_flow(V_test, theta_test, Y_bus_numerical)
hard_coded_time = time.time() - start_time
print(f"Performance comparison for {n_iterations} iterations:")
print(f"Symbolic-generated code time: {symbolic_time:.6f} seconds")
print(f"Hard-coded implementation time: {hard_coded_time:.6f} seconds")
print(f"Ratio (symbolic/hard-coded): {symbolic_time/hard_coded_time:.2f}")
# Example: Using symbolic-generated Jacobian
try:
J_numeric = jacobian_function(*all_vars_test)
print("\nJacobian evaluation successful")
print(f"Jacobian shape: {J_numeric.shape}")
# Simple verification: compare a few entries with numerical derivatives
print("\nVerification of Jacobian accuracy:")
# Calculate numerical derivatives for comparison
epsilon = 1e-6
# Choose a variable to perturb (e.g., first theta)
var_idx = n_bus # First theta variable (after V variables)
var_name = "theta_1"
# Make a copy of the test data
perturbed_vars = all_vars_test.copy()
# Perturb the variable
perturbed_vars[var_idx] += epsilon
# Calculate power with perturbed variable
perturbed_P = []
for i in range(min(2, len(power_functions['P']))):
perturbed_P.append(power_functions['P'][i](*perturbed_vars))
# Calculate numerical derivatives
numerical_derivs = []
for i in range(min(2, len(perturbed_P))):
original_P = power_functions['P'][i](*all_vars_test)
numerical_derivs.append((perturbed_P[i] - original_P) / epsilon)
# Compare with Jacobian entries
for i in range(min(2, len(numerical_derivs))):
symbolic_deriv = J_numeric[i, 0] if J_numeric.shape[1] > 0 else "N/A"
print(f"∂P_{i+1}/∂{var_name}:")
print(f" Numerical: {numerical_derivs[i]}")
print(f" Symbolic: {symbolic_deriv}")
except Exception as e:
print(f"Jacobian evaluation failed: {e}")
except Exception as e:
print(f"Performance comparison failed: {e}")
print("Using simplified example instead...")
# Simplified example with a single symbolic expression
x, y = sp.symbols('x y')
expr = sp.sin(x) * sp.cos(y) + x**2 * y
# Convert to numerical function
f_symbolic = lambdify([x, y], expr, 'numpy')
# Hard-coded equivalent
def f_hardcoded(x, y):
return np.sin(x) * np.cos(y) + x**2 * y
# Test values
x_val, y_val = 0.5, 1.0
# Time symbolic function
start_time = time.time()
for _ in range(10000):
result_symbolic = f_symbolic(x_val, y_val)
symbolic_time = time.time() - start_time
# Time hard-coded function
start_time = time.time()
for _ in range(10000):
result_hardcoded = f_hardcoded(x_val, y_val)
hardcoded_time = time.time() - start_time
print(f"Simplified example results:")
print(f" Symbolic result: {result_symbolic}")
print(f" Hard-coded result: {result_hardcoded}")
print(f" Symbolic time: {symbolic_time:.6f} s")
print(f" Hard-coded time: {hardcoded_time:.6f} s")
print(f" Ratio (symbolic/hard-coded): {symbolic_time/hardcoded_time:.2f}")
system = load_data('case3.xlsx')
# Generate power flow functions
try:
power_functions, jacobian_function = generate_power_flow_functions(system)
# Create test data
n_bus = system.bus.get_count()
V_test = np.ones(n_bus)
theta_test = np.zeros(n_bus)
all_vars_test = np.concatenate([V_test, theta_test])
# Evaluate power flow equations
print("Evaluating symbolic-generated power flow equations:")
for i in range(n_bus):
if i < len(power_functions['P']):
P_i = power_functions['P'][i](*all_vars_test)
print(f"P_{i+1} = {P_i}")
if i < len(power_functions['Q']):
Q_i = power_functions['Q'][i](*all_vars_test)
print(f"Q_{i+1} = {Q_i}")
# Test Jacobian function
try:
J_numeric = jacobian_function(*all_vars_test)
print(f"\nJacobian shape: {J_numeric.shape}")
if J_numeric.size > 0:
print("First Jacobian element:", J_numeric[0, 0])
except Exception as e:
print(f"Jacobian evaluation failed: {e}")
except Exception as e:
print(f"Function generation failed: {e}")
# Compare performance
performance_comparison(system)
Evaluating symbolic-generated power flow equations:
P_1 = 0
Q_1 = 0.5*b_1_2 + 0.5*b_3_1
P_2 = 50
Q_2 = 0.5*b_1_2 + 0.5*b_2_3
P_3 = -120
Q_3 = 0.5*b_2_3 + 0.5*b_3_1 - 50
Jacobian shape: (3, 3)
First Jacobian element: -1.0*x_1_2/(r_1_2**2 + x_1_2**2) - 1.0*x_2_3/(r_2_3**2 + x_2_3**2)
Performance comparison for 100 iterations:
Symbolic-generated code time: 0.129798 seconds
Hard-coded implementation time: 0.002911 seconds
Ratio (symbolic/hard-coded): 44.58
Jacobian evaluation successful
Jacobian shape: (3, 3)
Verification of Jacobian accuracy:
∂P_1/∂theta_1:
Numerical: -5.00044450291171e-7*r_1_2/(r_1_2**2 + x_1_2**2) - 5.00044450291171e-7*r_3_1/(r_3_1**2 + x_3_1**2) - 0.999999999999833*x_1_2/(r_1_2**2 + x_1_2**2) - 0.999999999999833*x_3_1/(r_3_1**2 + x_3_1**2)
Symbolic: -1.0*x_1_2/(r_1_2**2 + x_1_2**2) - 1.0*x_2_3/(r_2_3**2 + x_2_3**2)
∂P_2/∂theta_1:
Numerical: -5.00044450291171e-7*r_1_2/(r_1_2**2 + x_1_2**2) + 0.999999999999833*x_1_2/(r_1_2**2 + x_1_2**2)
Symbolic: 1.0*x_2_3/(r_2_3**2 + x_2_3**2)
Optimization Opportunities
When converting symbolic expressions to numerical code, there are several optimization opportunities:
Common subexpression elimination: Identify and reuse common intermediate results
Sparsity exploitation: Take advantage of zeros in the matrices to reduce computations
Parallelization: Independent calculations can be parallelized for improved performance
Vectorization: Use NumPy’s vectorized operations for efficiency
SymPy can help with some of these optimizations automatically, particularly common subexpression elimination.
With these techniques, we can enjoy the benefits of symbolic computation for analysis and insight while maintaining computational efficiency for practical applications.
Quizzes#
display_quiz(git_url + "sympy.json", colors=color_dict)