PINNs I: About stability and convergence

PINNs series - This is the first of a three-series blog posts about Physics-Informed Neural Networks (PINNs).

TL;DR. A generator will behave according to the swing equation, where its physical parameters and the magnitude of the disturbance will take it into a new state, possibly to a non-stable state too! We will solve this equation for a single machine system using PINNs $\rightarrow m \ddot{\delta}(t) + d \dot{\delta}(t) + B sin(\delta(t)) - P(t) = 0$

The future of the power system

Given the impact of the climate change crisis in our planet and the effort in the power sector to decarbonize the electricity supply, it is imperative to further increase the electricity produced by renewable resources.1 However, with a rising penetration of renewables and a decrease in conventional units, the system stability becomes more tricky since some important parameters of the system (inertia, explained below) will tend to decrease. Hence, the analysis of the stability of the power system becomes an important challenge for the future decarbonized grid.2

In general, the dynamic stability of the grid is analyzed by running many simulations and experiments based on ordinary differential equation (ODE) solvers. These solvers take some time to run and require tuning and special handling, this is where machine learning applications might come to the rescue to speed-up the calculations and simplify the process. However, keep in mind that the power system community has spent entire decades on developing these specialized tools! A great solution would take into account the knowledge we currently have of the power system and combine it with ML algorithms to further push the boundary of science 🤩

The swing equation

SMIB Figure 1. Single machine infinite bus system

First, let’s tap into the stability of power systems to study the physical dynamics of an electrical generator. The single machine infinite bus (SMIB) system helps us represent the behavior of a machine connected to the grid, as depicted in Figure 1.3 Here we have a generator $P_1$ connected to node 1, a transmission or distribution line represented by the susceptance $B_{12}$ between node 1 and node 2, and finally the power grid at node 2. Given some power demand and the difference in the voltages (driven by the voltage drop in the line) between both nodes, a power flows from the generator into the grid, represented by $P$. If we imagine this in the real world, $P_1$ could represent a wind turbine or gas turbine connected to some substation through a transmission line. This system can be represented by the swing equation, which directly originates from Newton’s second law for a rotational motion of a rigid body:

\[T_a(t) = J \alpha_m(t)\] \[T_a(t) = T_m(t) - T_e(t)\]

where

$T \rightarrow$ (a)ccelerating, (m)echanical and (e)lectrical torque [$\text{N-m}$]

$J \rightarrow$ moment of inertia or rotational inertia [$\text{kg}-m^2$]

$\alpha \rightarrow$ rotor angular acceleration [$\text{rad}/s^2$]

This equation describes how, given a difference between a mechanical torque (a force given by some driver like a turbine connected to the generator) and an electrical torque (a force requested by the power grid), the generator will respond with certain acceleration $\alpha_m$. The moment of inertia or rotational inertia can be interpreted as the proportional of the torque that needs to be given to the system to start accelerating (or decelerating). In simpler terms, it basically tells us that if we apply torque to a rotational body, then it will start to accelerate if we give it enough force (or the other way around, if the body is already rotating the moment of inertia opposes to change, so you need to give some more torque in the opposite direction to start decelerating). Notice that all terms depend on time except the inertia, which is a property of the body!

Further, the angular acceleration of the rotor can be decomposed into its derivatives:

\[\alpha_m(t) = \frac{dw_m(t)}{dt} = \frac{d^2\theta_m(t)}{dt^2}\] \[w_m(t) = \frac{d\theta_m(t)}{dt}\]

where

$w \rightarrow$ rotor angular velocity [$\text{rad}/s$]

$\theta \rightarrow$ rotor angle [$\text{rad}$]

In steady-state, $T_m = T_e$ which makes $T_a = 0$, i.e. the rotor acceleration is zero and there’s a constant rotor velocity. This means that in a stable state there’s no acceleration, just a constant rotor velocity. Imagine that if two people are pulling a string, one in each side, if they both exert the exact same force they will be in an equilibrium state, not accelerating at all! However, in a rotational body the system keeps a constant angular velocity, which is basically the velocity at which the rotor rotates (that’s why it is called rotor!).

Now, if we define the synchronous speed as the reference axis, which is the speed at which the magnetic field rotates to keep some designed frequency (50 Hz in Europe, 60 Hz in North America). And if we magically had the ability to measure the rotor angular position at every time step with respect to this reference axis, then:

\[\theta_m(t) = w_{syn}(t) + \delta_m(t)\]

where

$w_{syn} \rightarrow$ synchronous angular velocity of rotor [$\text{rad}/s$]

$\delta_m \rightarrow$ rotor angular position w.r.t. synchronous reference [$\text{rad}$]

Then using these definitions in the swing equation:

\[T_m(t) - T_e(t) = J\frac{d^2\theta_m(t)}{dt^2}\]

Measuring only from the reference axis:

\[T_m(t) - T_e(t) = J \frac{d^2\delta_m(t)}{dt^2}\]

Multiplying both sides by $w_m$ (to obtain power instead of torque) and dividing by $\text{S}_{rated}$ (to obtain power in per unit):

\[p_{m,}(t) - p_{e}(t) = \frac{Jw_m(t)}{\text{S}_{rated}}\frac{d^2\delta(t)}{dt^2}\]

Adding a damping term $d$, which is similar to the rotational inertia but depends on the angular velocity:

\[p_{m}(t) - p_{e} = \frac{Jw_m(t)}{\text{S}_{rated}}\frac{d^2\delta(t)}{dt^2} + d\frac{d\delta(t)}{dt}\]

Defining an inertia constant as $m$:

\[p_{m}(t) - p_{e}(t) = m\frac{d^2\delta(t)}{dt^2} + d\frac{d\delta(t)}{dt}\]

Finally, using Newton’s dot notation for derivatives:

\[p_{m}(t) - p_{e}(t) = m\ddot{\delta}(t) + d\dot{\delta}(t)\]

Notice. This formulation ignores transmission losses and bus voltage deviations!

ODE Formulation

Now that we defined the basis of the swing equation, we can represent the SMIB system as an ordinary differential equation (ODE):

\[m_k \ddot{\delta} + d_k \dot{\delta} + B_{kj} V_k V_j sin(\delta) - P_k = 0\]

Then, we can express this second-order differential equation into a system of ODEs by defining the state-vector:

\[\begin{bmatrix}x_1 \\ x_2\end{bmatrix} = \begin{bmatrix} \delta\\ \dot{\delta} \end{bmatrix}\]

where

\[\frac{d}{dt}\begin{bmatrix} \delta \\ \dot{\delta} \end{bmatrix} = \begin{bmatrix} \dot{\delta} \\ \ddot{\delta} \end{bmatrix}\] \[\begin{bmatrix} \dot{\delta} \\ \ddot{\delta} \end{bmatrix} = \begin{bmatrix} \dot{\delta} \\ \frac{P - d \dot{\delta} - B\ \text{sin}\left( \delta_i - \delta_j \right)}{m} \end{bmatrix}\]

and, since $\delta_j$ is always zero in our case:

\[\begin{bmatrix} \dot{\delta} \\ \ddot{\delta} \end{bmatrix} = \begin{bmatrix} \dot{\delta} \\ \frac{P - d \dot{\delta} - B\ \text{sin}\left( \delta \right)}{m} \end{bmatrix}\]

In code, it would look like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def swing_equation(t, x, m, d, B, P):
    r""" Swing equation expressed as a system of ODEs.
    
    Given by the formula:
    
    .. math::
        m_i \ddot{\delta} + d_i \dot{\delta} + B_{ij} V_i V_j sin(\delta) - P_i = 0
    
    Args:
        x: State vector [delta_i, omega_i]
        m: Inertia of the machine
        d: Damping coefficient
        B: Bus susceptance matrix
        P: Power injection or retrieval
    
    Returns:
        x_new: Updated state vector [delta_i, omega_i]
    """
    # Split the state variable into delta and omega
    state_delta = x[0]
    state_omega = x[1]

    # Computing the non-linear term in the swing equation sum_j (B_ij sin(delta_i - delta_j))
    delta_i = state_delta
    delta_j = np.zeros_like(delta_i)
    delta_ij = np.sin(delta_i - delta_j)
    connectivity_vector = np.sum(np.multiply(B, delta_ij))

    # Preallocate memory
    state_delta_new = np.zeros_like(state_delta)
    state_omega_new = np.zeros_like(state_omega)
    
    # Update states
    state_delta_new = state_omega
    state_omega_new = 1 / m * (P - d * state_omega - connectivity_vector)

    return np.array([state_delta_new, state_omega_new])

And to solve it we can use the integrate function from scipy, as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def solve_swing_eq(params):
    """Solves the swing equation with pre-defined parameters.
    
    Args:
        params: Dictionary containing relevant parameters.
    
    Returns:
        states: Solution of the swing equation (angle)
    """

    # Unpack parameters
    m = params['inertia']
    d = params['damping']
    B = params['susceptance']
    P = params['power']
    delta_0 = np.array([params['delta_0']])
    omega_0 = np.array([params['omega_0']])
    t_span = params['t_span'] # (t_min, t_max)
    n_data = params['n_data']
    
    # Define analysis time period
    t_eval = np.linspace(t_span[0], t_span[1], n_data)
    
    # Pack initial state
    initial_states = np.concatenate([delta_0, omega_0])
    
    # Solve ODE
    ode_solution = integrate.solve_ivp(swing_equation,
                                       t_span=t_span,
                                       y0=initial_states,
                                       args=[m, d, B, P],
                                       t_eval=t_eval)
    
    # Extract solution only
    states = ode_solution.y
    
    return states

These two code chunks work for a time span, but for a specific power disturbance $P$. We would need to solve it for as many power disturbances cases as we want.

About stability and convergence

By solving the previous ODE over time, we can obtain the exact behavior of the machine given some power disturbance. Let’s analyze the case when $P = 0.01$, $P = 0.1$, $P = 0.2$, and $P = 0.3$.

Rotor angle Figure 2. Rotor angle vs. Time

On the rotor angle, in the first two cases it seems that the angle increases with the disturbance to reach an apex and a nadir in the first cycle, then stabilizes. Notice the y-axis, with a very small disturbance the angle doesn’t change much at 0.05 radians, while with a disturbance of 0.1 the angle increases by 0.5 radians. Now, the case with $P=0.3$ is quite interesting, in this case the ODE enters a state of unstability and diverges, i.e. increases the angle and doesn’t converge to a certain value.

Rotor speed Figure 3. Rotor speed vs. Time

We can see the same behavior on the rotational speed of the rotor, in the first two cases the value stabilizes but on the third case it follows a constant cyclical pattern. Remember that the derivative of the speed is the angle, that’s why the previous plot increased with small bumps, these are shown here as a cyclical pattern.

Now, let’s take a look at the phase plots, same cases:

Phase plot Figure 4. Phase plot - Rotor speed vs. Rotor angle

Here we can see the derivatives as an arrow, given a specific point (at certain $\delta$ and $\omega$). In the first two cases, the gradient takes the system into a stable point, starting at 0. However, in the last case, the system follows the gradients and stays in this cyclical pattern we saw in the previous two plots.

Final remarks

In this first post, the swing equation case was proposed with an initial investigation on the convergence of the system. This will help us interpret the behavior of the PINNs given this stable/non-stable case of the equation.

References

1: D. Rolnick. et. al., Tackling Climate Change with Machine Learning. Paper

2: P. Denholm, T. Mai, R. Wallace Kenyon, B. Kroposki, M. O’Malley. Inertia and the Power Grid: A Guide Without the Spin. National Renewable Energy Laboratory (NREL). Report

3: G. S. Misyris, A. Venzke, S. Chatzivasileiadis, Physics-Informed Neural Networks for Power Systems, accepted at IEEE PES General Meeting 2020, Montreal, Canada. Paper

Updated:

Leave a comment