From 00273eddc78b0a2e2272363fa089d21cfd83ce83 Mon Sep 17 00:00:00 2001
From: Marko Hinkkanen <76600872+mhinkkan@users.noreply.github.com>
Date: Sun, 8 Sep 2024 13:23:05 +0300
Subject: [PATCH] Add grid-forming control documentation, minor refactoring
 (#167)

* Add documentation about grid-forming control plus minor refactoring.

* Trimming white space
---
 docs/source/control/control_system.rst        | 30 +++---
 docs/source/control/drive/current_ctrl.rst    | 46 +++++-----
 docs/source/control/drive/observers.rst       | 36 ++++----
 docs/source/control/grid/current_ctrl.rst     | 25 +++--
 docs/source/control/grid/dc_voltage_ctrl.rst  | 16 ++--
 docs/source/control/grid/grid-forming.rst     | 92 +++++++++++++++++++
 docs/source/control/grid/index.rst            |  5 +-
 .../grid/{synchronization.rst => pll.rst}     | 57 +++++++-----
 examples/drive/flux_vector/README.txt         |  2 +-
 .../drive/vhz/plot_vhz_ctrl_6step_im_2kw.py   |  8 +-
 examples/drive/vhz/plot_vhz_ctrl_im_2kw.py    |  2 +-
 .../grid/grid_following/plot_gfl_10kva.py     |  4 +-
 .../grid/grid_forming/plot_gfm_obs_13kva.py   |  6 +-
 .../grid/grid_forming/plot_gfm_rfpsc_13kva.py |  9 +-
 motulator/drive/control/im/_vhz.py            |  6 +-
 motulator/drive/control/sm/_flux_vector.py    | 30 +++---
 motulator/grid/control/_common.py             | 46 +++++-----
 motulator/grid/control/_grid_following.py     |  3 +
 motulator/grid/control/_observer_gfm.py       | 81 ++++++++--------
 .../grid/control/_power_synchronization.py    | 21 +++--
 motulator/grid/utils/_plots.py                | 13 +--
 pyproject.toml                                |  1 +
 22 files changed, 317 insertions(+), 222 deletions(-)
 create mode 100644 docs/source/control/grid/grid-forming.rst
 rename docs/source/control/grid/{synchronization.rst => pll.rst} (78%)

diff --git a/docs/source/control/control_system.rst b/docs/source/control/control_system.rst
index 4b3c4c3b..a725ccbb 100644
--- a/docs/source/control/control_system.rst
+++ b/docs/source/control/control_system.rst
@@ -4,11 +4,11 @@ Common
 Data Flow
 ---------
 
-The two figures below show the structure and data flow in a typical simulation model as well as an example of the internal structure of a discrete-time control system. The text in italics refers to the default object names used in the software. In discrete-time control systems, the signals are collected into the following :class:`types.SimpleNamespace` objects: 
+The two figures below show the structure and data flow in a typical simulation model as well as an example of the internal structure of a discrete-time control system. The text in italics refers to the default object names used in the software. In discrete-time control systems, the signals are collected into the following :class:`types.SimpleNamespace` objects:
 
-- The object `fbk` contains feedback signals for the controllers. These signals can be measured quantities (such as the measured DC-bus voltage `fbk.u_dc`) as well as estimated quantities (such as the estimated stator flux linkage `fbk.psi_s`). 
+- The object `fbk` contains feedback signals for the controllers. These signals can be measured quantities (such as the measured DC-bus voltage `fbk.u_dc`) as well as estimated quantities (such as the estimated stator flux linkage `fbk.psi_s`).
 
-- The object `ref` contains the reference signals, both externally provided (such as `ref.w_m` for the angular speed reference) and internally generated (such as `ref.d_c_abc` for the duty ratios). 
+- The object `ref` contains the reference signals, both externally provided (such as `ref.w_m` for the angular speed reference) and internally generated (such as `ref.d_c_abc` for the duty ratios).
 
 These objects `fbk` and `ref` may propagate through several blocks (implemented as classes in the software), which may add new signals or modify the existing one. At the end of their propagation chain, both objects `fbk` and `ref` are saved to the lists. Therefore, the most relevant internal signals of the control system are directly accessible after the simulation. Furthermore, the states and inputs of the continuous-time system model are also saved. In the post-processing stage, the saved data is converted to the NumPy arrays.
 
@@ -16,22 +16,22 @@ These objects `fbk` and `ref` may propagate through several blocks (implemented
    :width: 100%
    :align: center
    :alt: Block diagram of the control system.
-   
-   Block diagram illustrating the structure and data flow in a typical simulation model. The lower part of the figure illustrates how the data is saved. The post-processing is done after the simulation. The internal structure of a typical control system is exemplified in the figure below. 
+
+   Block diagram illustrating the structure and data flow in a typical simulation model. The lower part of the figure illustrates how the data is saved. The post-processing is done after the simulation. The internal structure of a typical control system is exemplified in the figure below.
 
 .. figure:: figs/discrete_control_system.svg
    :width: 100%
    :align: center
    :alt: Block diagram of the control system.
-   
-   Block diagram exemplifying the internal structure of a typical cascade control system. The object `ref` at the control system output should contain the sampling period `T_s` and the converter duty ratios `d_c_abc` for the carrier comparison. The observer does not necessarily exist in all control systems (or it can be replaced with, e.g., a phase-locked loop). 
+
+   Block diagram exemplifying the internal structure of a typical cascade control system. The object `ref` at the control system output should contain the sampling period `T_s` and the converter duty ratios `d_c_abc` for the carrier comparison. The observer does not necessarily exist in all control systems (or it can be replaced with, e.g., a phase-locked loop).
 
 Main Control Loop
 -----------------
 
 A template for the main control loop is available in the base class for discrete-time control systems in :class:`motulator.common.control.ControlSystem`. The main control loop in this template has the following steps:
 
-1. Get the feedback signals `fbk` for the controllers from the outputs of the continuous-time system model `mdl`. This step may contain first getting the measurements and then optionally computing the observer outputs (or otherwise manipulating the measured signals).  
+1. Get the feedback signals `fbk` for the controllers from the outputs of the continuous-time system model `mdl`. This step may contain first getting the measurements and then optionally computing the observer outputs (or otherwise manipulating the measured signals).
 2. Get the reference signals `ref` and compute the controller outputs based on the feedback signals `fbk`. Cascade control systems may contain multiple controllers, where the output of the outer controller is the reference signal for the inner controller.
 3. Update the states of the control system for the next sampling instant.
 4. Save the feedback signals `fbk` and the reference signals `ref` so they can be accessed after the simulation.
@@ -54,14 +54,14 @@ The figure below shows a 2DOF PI controller with an optional feedforward term. I
    u &= k_\mathrm{t}r - k_\mathrm{p}y + u_\mathrm{i} + u_\mathrm{ff}
    :label: 2dof_pi
 
-where :math:`r` is the reference signal, :math:`y` is the measured (or estimated) feedback signal, :math:`u_\mathrm{i}` is the the integral state, and :math:`u_\mathrm{ff}` is the optional feedforward signal. Furthermore, :math:`k_\mathrm{t}` is the reference-feedforward gain, :math:`k_\mathrm{p}` is the proportional gain, and :math:`k_\mathrm{i}` is the integral gain. Setting :math:`k_\mathrm{t} = k_\mathrm{p}` and :math:`u_\mathrm{ff} = 0` results in the standard PI controller. This 2DOF PI controller can also be understood as a state-feedback controller with integral action and reference feedforward [#Fra1997]_. 
+where :math:`r` is the reference signal, :math:`y` is the measured (or estimated) feedback signal, :math:`u_\mathrm{i}` is the the integral state, and :math:`u_\mathrm{ff}` is the optional feedforward signal. Furthermore, :math:`k_\mathrm{t}` is the reference-feedforward gain, :math:`k_\mathrm{p}` is the proportional gain, and :math:`k_\mathrm{i}` is the integral gain. Setting :math:`k_\mathrm{t} = k_\mathrm{p}` and :math:`u_\mathrm{ff} = 0` results in the standard PI controller. This 2DOF PI controller can also be understood as a state-feedback controller with integral action and reference feedforward [#Fra1997]_.
 
 .. figure:: figs/2dof_pi.svg
    :width: 100%
    :align: center
    :alt: 2DOF PI controller.
 
-   2DOF PI controller with an optional feedforward term. The operator :math:`1/s` refers to integration. A discrete-time variant of this controller with the integrator anti-windup is implemented in the :class:`motulator.common.control.PIController` class.  
+   2DOF PI controller with an optional feedforward term. The operator :math:`1/s` refers to integration. A discrete-time variant of this controller with the integrator anti-windup is implemented in the :class:`motulator.common.control.PIController` class.
 
 Disturbance-Observer Structure
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -71,7 +71,7 @@ The controller :eq:`2dof_pi` can be equally represented using the disturbance-ob
 .. math::
 	\frac{\mathrm{d} u_\mathrm{i}}{\mathrm{d} t} &= \alpha_\mathrm{i}\left(u - \hat v\right) \\
    \hat v &= u_\mathrm{i} - (k_\mathrm{p} - k_\mathrm{t})y + u_\mathrm{ff} \\
-   u &= k_\mathrm{t}\left(r - y\right) + \hat v 
+   u &= k_\mathrm{t}\left(r - y\right) + \hat v
    :label: 2dof_pi_disturbance_observer
 
 where :math:`\alpha_\mathrm{i} = k_\mathrm{i}/k_\mathrm{t}` is the redefined integral gain and :math:`\hat v` is the input-equivalent disturbance estimate. This structure is convenient to prevent the integral windup that originates from the actuator saturation [#Fra1997]_. The actuator output is limited in practice due to physical constraints. Consequently, the realized actuator output is
@@ -85,7 +85,7 @@ where :math:`\mathrm{sat}(\cdot)` is the saturation function. If this saturation
 	\frac{\mathrm{d} u_\mathrm{i}}{\mathrm{d} t} = \alpha_\mathrm{i}\left(\bar{u} - \hat v \right)
    :label: anti_windup
 
-The other parts of the above controller are not affected by the saturation. 
+The other parts of the above controller are not affected by the saturation.
 
 Discrete-Time Algorithm
 ^^^^^^^^^^^^^^^^^^^^^^^
@@ -99,7 +99,7 @@ The discrete-time variant of the controller :eq:`2dof_pi_disturbance_observer` w
    \bar{u}(k) &= \mathrm{sat}[u(k)]
    :label: discrete_2dof_pi
 
-where :math:`T_\mathrm{s}` is the sampling period and :math:`k` is the discrete-time index. This algorithm corresponds to the actual implementation in the :class:`motulator.common.control.PIController` class. 
+where :math:`T_\mathrm{s}` is the sampling period and :math:`k` is the discrete-time index. This algorithm corresponds to the actual implementation in the :class:`motulator.common.control.PIController` class.
 
 .. _complex-vector-2dof-pi-controller:
 
@@ -113,7 +113,7 @@ As shown in the figure below, the 2DOF PI controller presented above can be exte
     \boldsymbol{u} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{r} - \boldsymbol{k}_\mathrm{p}\boldsymbol{y} + \boldsymbol{u}_\mathrm{i} + \boldsymbol{u}_\mathrm{ff}
     :label: complex_vector_2dof_pi
 
-where :math:`\boldsymbol{u}` is the output of the controller, :math:`\boldsymbol{r}` is the reference signal, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\boldsymbol{u}_\mathrm{ff}` is the optional feedforward signal. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain. 
+where :math:`\boldsymbol{u}` is the output of the controller, :math:`\boldsymbol{r}` is the reference signal, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\boldsymbol{u}_\mathrm{ff}` is the optional feedforward signal. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain.
 
 .. figure:: figs/complex_vector_2dof_pi.svg
    :width: 100%
@@ -131,5 +131,3 @@ The discrete-time implementation of :eq:`complex_vector_2dof_pi` with the anti-w
 .. [#Fra1997] Franklin, Powell, Workman, "Digital Control of Dynamic Systems," 3rd ed., Menlo Park, CA: Addison-Wesley, 1997
 
 .. [#Bri1999] Briz del Blanco, Degner, Lorenz, “Dynamic analysis of current regulators for AC motors using complex vectors,” IEEE Trans. Ind. Appl., 1999, https://doi.org/10.1109/28.806058
-
-
diff --git a/docs/source/control/drive/current_ctrl.rst b/docs/source/control/drive/current_ctrl.rst
index 4a6dd676..98ba18a1 100644
--- a/docs/source/control/drive/current_ctrl.rst
+++ b/docs/source/control/drive/current_ctrl.rst
@@ -1,7 +1,7 @@
 Current Control
 ===============
 
-Synchronous-frame two-degrees-of-freedom (2DOF) proportional-integral (PI) current control is commonly used in three-phase AC machine drives [#Har1998]_, [#Bri1999]_, [#Awa2019]_, [#Hin2024]_. This control structure allows to compensate for the cross-coupling originating from rotating coordinates as well as to improve disturbance rejection. 
+Synchronous-frame two-degrees-of-freedom (2DOF) proportional-integral (PI) current control is commonly used in three-phase AC machine drives [#Har1998]_, [#Bri1999]_, [#Awa2019]_, [#Hin2024]_. This control structure allows to compensate for the cross-coupling originating from rotating coordinates as well as to improve disturbance rejection.
 
 A 2DOF PI current controller for induction machines is available in the :class:`motulator.drive.control.im.CurrentController` class and for synchronous machines in the :class:`motulator.drive.control.sm.CurrentController` class, both of which inherit from the :class:`motulator.common.control.ComplexPIController` class. In the following, current control of induction machines is first considered in detail. Then, the same principles are applied to synchronous machines. Complex space vectors are used to represent three-phase quantities.
 
@@ -15,7 +15,7 @@ The inverse-Γ model of an induction machine is considered (see :doc:`/model/dri
 
 .. math::
     L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t} &= \boldsymbol{u}_\mathrm{s} - (R_\sigma + \mathrm{j} \omega_\mathrm{s}L_\sigma)\boldsymbol{i}_\mathrm{s} - \underbrace{\left(\mathrm{j}\omega_\mathrm{m} - \frac{R_\mathrm{R}}{L_\mathrm{M}}\right)\boldsymbol{\psi}_\mathrm{R}}_{\text{back-emf } \boldsymbol{e}_\mathrm{s}} \\
-	\frac{\mathrm{d} \boldsymbol{\psi}_\mathrm{R}}{\mathrm{d} t} &= R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\frac{R_\mathrm{R}}{L_\mathrm{M}} + \mathrm{j}\omega_\mathrm{r} \right)\boldsymbol{\psi}_\mathrm{R} 
+	\frac{\mathrm{d} \boldsymbol{\psi}_\mathrm{R}}{\mathrm{d} t} &= R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\frac{R_\mathrm{R}}{L_\mathrm{M}} + \mathrm{j}\omega_\mathrm{r} \right)\boldsymbol{\psi}_\mathrm{R}
 
 where :math:`R_\sigma = R_\mathrm{s} + R_\mathrm{R}` is the total resistance and :math:`\omega_\mathrm{r} = \omega_\mathrm{s} - \omega_\mathrm{m}` is the slip angular frequency. The rotor flux linkage :math:`\boldsymbol{\psi}_\mathrm{R}` and the rotor speed :math:`\omega_\mathrm{m}` change slowly as compared to the stator current. Consequently, the back-emf :math:`\boldsymbol{e}_\mathrm{s}` can be considered as a quasi-constant input (load) disturbance for the current controller, and it suffices to consider the stator current dynamics
 
@@ -30,15 +30,15 @@ The design of synchronous-frame 2DOF PI current control is considered in the con
 
 .. math::
 	\frac{\mathrm{d} \boldsymbol{u}_\mathrm{i}}{\mathrm{d} t} &= (\boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s}\boldsymbol{k}_\mathrm{t} )\left(\boldsymbol{i}_\mathrm{s,ref} - \boldsymbol{i}_\mathrm{s}\right) \\
-    \boldsymbol{u}_\mathrm{s,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{i}_\mathrm{s,ref} - \boldsymbol{k}_\mathrm{p}\boldsymbol{i}_\mathrm{s} + \boldsymbol{u}_\mathrm{i} 
+    \boldsymbol{u}_\mathrm{s,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{i}_\mathrm{s,ref} - \boldsymbol{k}_\mathrm{p}\boldsymbol{i}_\mathrm{s} + \boldsymbol{u}_\mathrm{i}
     :label: cc
 
-where :math:`\boldsymbol{u}_\mathrm{s,ref}` is the output of the controller, i.e., the stator voltage reference, :math:`\boldsymbol{i}_\mathrm{s,ref}` is the stator current reference, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\omega_\mathrm{s}` is the angular speed of the coordinate system. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain. 
+where :math:`\boldsymbol{u}_\mathrm{s,ref}` is the output of the controller, i.e., the stator voltage reference, :math:`\boldsymbol{i}_\mathrm{s,ref}` is the stator current reference, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\omega_\mathrm{s}` is the angular speed of the coordinate system. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain.
 
 .. note::
-   The gain definitions used in :eq:`cc` differ from [#Hin2024]_, where a more general controller structure is considered. 
+   The gain definitions used in :eq:`cc` differ from [#Hin2024]_, where a more general controller structure is considered.
 
-Closed-Loop System 
+Closed-Loop System
 ^^^^^^^^^^^^^^^^^^
 
 Here, ideal voltage production is assumed, :math:`\boldsymbol{u}_\mathrm{s} = \boldsymbol{u}_\mathrm{s,ref}`. Using :eq:`im_current` and :eq:`cc`, the closed-loop system in the Laplace domain becomes
@@ -49,13 +49,13 @@ Here, ideal voltage production is assumed, :math:`\boldsymbol{u}_\mathrm{s} = \b
 The disturbance rejection depends on the closed-loop admittance
 
 .. math::
-    \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L_\sigma s^2 + (R_\sigma + \mathrm{j}\omega_\mathrm{s} L_\sigma + \boldsymbol{k}_\mathrm{p}) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s} \boldsymbol{k}_\mathrm{t}} 
+    \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L_\sigma s^2 + (R_\sigma + \mathrm{j}\omega_\mathrm{s} L_\sigma + \boldsymbol{k}_\mathrm{p}) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s} \boldsymbol{k}_\mathrm{t}}
     :label: Yc
 
 The closed-loop poles can be arbitrarily placed by means of the gains. The reference-tracking transfer function is
 
 .. math::
-	\boldsymbol{G}_\mathrm{c}(s) = \frac{(s + \mathrm{j}\omega_\mathrm{s}) \boldsymbol{k}_\mathrm{t} + \boldsymbol{k}_\mathrm{i} }{L_\sigma s^2 + (R_\sigma + \mathrm{j}\omega_\mathrm{s} L_\sigma + \boldsymbol{k}_\mathrm{p}) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s} \boldsymbol{k}_\mathrm{t}}     
+	\boldsymbol{G}_\mathrm{c}(s) = \frac{(s + \mathrm{j}\omega_\mathrm{s}) \boldsymbol{k}_\mathrm{t} + \boldsymbol{k}_\mathrm{i} }{L_\sigma s^2 + (R_\sigma + \mathrm{j}\omega_\mathrm{s} L_\sigma + \boldsymbol{k}_\mathrm{p}) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s} \boldsymbol{k}_\mathrm{t}}
     :label: Gc
 
 whose zero can be placed by means of the reference-feedforward gain :math:`\boldsymbol{k}_\mathrm{t}`.
@@ -65,7 +65,7 @@ Gain Selection
 
 Consider the gains
 
-.. math::                
+.. math::
     \boldsymbol{k}_\mathrm{p} = 2\alpha_\mathrm{c} \hat L_\sigma - \hat R_\sigma \qquad\qquad
     \boldsymbol{k}_\mathrm{i} = \alpha_\mathrm{c}^2\hat L_\sigma  \qquad \qquad
     \boldsymbol{k}_\mathrm{t} = \alpha_\mathrm{c} \hat L_\sigma
@@ -76,9 +76,9 @@ where :math:`\hat R_\sigma = 0` can be used in practice due to its minor effect
 .. math::
     \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L_\sigma (s + \alpha_\mathrm{c})(s + \alpha_\mathrm{c} + \mathrm{j}\omega_\mathrm{s})}
     \qquad\qquad
-    \boldsymbol{G}_\mathrm{c}(s) = \frac{\alpha_\mathrm{c}}{s + \alpha_\mathrm{c}} 
+    \boldsymbol{G}_\mathrm{c}(s) = \frac{\alpha_\mathrm{c}}{s + \alpha_\mathrm{c}}
 
-It can be seen that this design results in the first-order reference-tracking dynamics. Furthermore, one pole is placed at the real axis at :math:`s=-\alpha_\mathrm{c}`, while another pole moves with the angular frequency of the coordinate system, :math:`s= -\alpha_\mathrm{c} - \mathrm{j}\omega_\mathrm{s}`. The complex-vector design tends to be slightly more robust to parameter errors than the IMC design since the other closed-loop pole approximately corresponds to the open-loop pole.  
+It can be seen that this design results in the first-order reference-tracking dynamics. Furthermore, one pole is placed at the real axis at :math:`s=-\alpha_\mathrm{c}`, while another pole moves with the angular frequency of the coordinate system, :math:`s= -\alpha_\mathrm{c} - \mathrm{j}\omega_\mathrm{s}`. The complex-vector design tends to be slightly more robust to parameter errors than the IMC design since the other closed-loop pole approximately corresponds to the open-loop pole.
 
 This gain selection is used in the :class:`motulator.drive.control.im.CurrentController` class. The stator voltage is limited in practice due to the limited DC-bus voltage of the converter. Consequently, the realized (limited) voltage reference is
 
@@ -93,13 +93,13 @@ Synchronous Machines
 System Model
 ^^^^^^^^^^^^
 
-Consider the synchronous machine model in rotor coordinates, rotating at :math:`\omega_\mathrm{m}` and aligned along the d-axis of the rotor, 
+Consider the synchronous machine model in rotor coordinates, rotating at :math:`\omega_\mathrm{m}` and aligned along the d-axis of the rotor,
 
 .. math::
     \frac{\mathrm{d}\boldsymbol{\psi}_\mathrm{s}}{\mathrm{d} t} &= \boldsymbol{u}_\mathrm{s} - R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - \mathrm{j}\omega_\mathrm{m}\boldsymbol{\psi}_\mathrm{s} \\
- 	\boldsymbol{\psi}_\mathrm{s} &= L_\mathrm{d}\mathrm{Re}\{\boldsymbol{i}_\mathrm{s}\} + \mathrm{j}L_\mathrm{q}\mathrm{Im}\{\boldsymbol{i}_\mathrm{s}\} + \psi_\mathrm{f} 
+ 	\boldsymbol{\psi}_\mathrm{s} &= L_\mathrm{d}\mathrm{Re}\{\boldsymbol{i}_\mathrm{s}\} + \mathrm{j}L_\mathrm{q}\mathrm{Im}\{\boldsymbol{i}_\mathrm{s}\} + \psi_\mathrm{f}
     :label: sm_model
- 
+
 where linear magnetics are assumed for simplicity. For further details, see :doc:`/model/drive/machines`.
 
 2DOF PI Controller
@@ -109,28 +109,28 @@ An internal change of the state variable from the stator current to the stator f
 
 .. math::
     \boldsymbol{\psi}_\mathrm{s,ref} &= \hat{L}_\mathrm{d}\mathrm{Re}\{\boldsymbol{i}_\mathrm{s,ref}\} + \mathrm{j} \hat{L}_\mathrm{q}\mathrm{Im}\{\boldsymbol{i}_\mathrm{s,ref}\} \\
-    \hat{\boldsymbol{\psi}_\mathrm{s}} &= \hat{L}_\mathrm{d}\mathrm{Re}\{\boldsymbol{i}_\mathrm{s}\} + \mathrm{j} \hat{L}_\mathrm{q}\mathrm{Im}\{\boldsymbol{i}_\mathrm{s}\} 
+    \hat{\boldsymbol{\psi}_\mathrm{s}} &= \hat{L}_\mathrm{d}\mathrm{Re}\{\boldsymbol{i}_\mathrm{s}\} + \mathrm{j} \hat{L}_\mathrm{q}\mathrm{Im}\{\boldsymbol{i}_\mathrm{s}\}
     :label: flux_mapping_sm
 
-This choice of using the flux linkage as the internal state has some advantages: the gain expressions become simpler; the magnetic saturation would be more convenient to take into account; and the same control structure can be used for salient and nonsalient machines. 
+This choice of using the flux linkage as the internal state has some advantages: the gain expressions become simpler; the magnetic saturation would be more convenient to take into account; and the same control structure can be used for salient and nonsalient machines.
 
-Here, the complex vector design is considered. Hence, the controller :eq:`cc` can be rewritten as 
+Here, the complex vector design is considered. Hence, the controller :eq:`cc` can be rewritten as
 
 .. math::
 	\frac{\mathrm{d} \boldsymbol{u}_\mathrm{i}}{\mathrm{d} t} &= (\boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{s}\boldsymbol{k}_\mathrm{t} )\left(\boldsymbol{\psi}_\mathrm{s,ref} - \hat{\boldsymbol{\psi}}_\mathrm{s}\right) \\
-    \boldsymbol{u}_\mathrm{s,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{\psi}_\mathrm{s,ref} - \boldsymbol{k}_\mathrm{p}\hat{\boldsymbol{\psi}}_\mathrm{s} + \boldsymbol{u}_\mathrm{i} 
+    \boldsymbol{u}_\mathrm{s,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{\psi}_\mathrm{s,ref} - \boldsymbol{k}_\mathrm{p}\hat{\boldsymbol{\psi}}_\mathrm{s} + \boldsymbol{u}_\mathrm{i}
     :label: cc_flux
 
-where the angular speed of the coordinate system equals typically the measured rotor speed, :math:`\omega_\mathrm{s} = \omega_\mathrm{m}`, or the estimated rotor speed :math:`\omega_\mathrm{s} = \hat{\omega}_\mathrm{m}`. If the magnetic saturation is not considered, this flux-linkage-based current controller is equivalent to a regular 2DOF PI current controller (even if inductance estimates are inaccurate). Notice that :math:`\boldsymbol{i}_\mathrm{s,ref} = \boldsymbol{i}_\mathrm{s}` holds in the steady state even with inductance estimate inaccuracies, since the same inductances are used to map both the reference current and the actual current to the corresponding flux linkages. 
+where the angular speed of the coordinate system equals typically the measured rotor speed, :math:`\omega_\mathrm{s} = \omega_\mathrm{m}`, or the estimated rotor speed :math:`\omega_\mathrm{s} = \hat{\omega}_\mathrm{m}`. If the magnetic saturation is not considered, this flux-linkage-based current controller is equivalent to a regular 2DOF PI current controller (even if inductance estimates are inaccurate). Notice that :math:`\boldsymbol{i}_\mathrm{s,ref} = \boldsymbol{i}_\mathrm{s}` holds in the steady state even with inductance estimate inaccuracies, since the same inductances are used to map both the reference current and the actual current to the corresponding flux linkages.
 
 The gain selection analogous to :eq:`complex_vector_gains` becomes
 
-.. math::                
+.. math::
     \boldsymbol{k}_\mathrm{p} = 2\alpha_\mathrm{c} - \hat R_\mathrm{s} \qquad\qquad
     \boldsymbol{k}_\mathrm{i} = \alpha_\mathrm{c}^2 \qquad \qquad
-    \boldsymbol{k}_\mathrm{t} = \alpha_\mathrm{c} 
+    \boldsymbol{k}_\mathrm{t} = \alpha_\mathrm{c}
 
-where :math:`\hat R_\mathrm{s} = 0` can be used in practice. Assume accurate parameter estimates and perfect alignment of the controller coordinate system with the rotor coordinate system. Then, using :eq:`sm_model`, :eq:`flux_mapping_sm`, and :eq:`cc_flux`, the closed-loop system can be shown to be analogous to the induction machine case. This control design corresponds to the implementation in the :class:`motulator.drive.control.sm.CurrentController` class. 
+where :math:`\hat R_\mathrm{s} = 0` can be used in practice. Assume accurate parameter estimates and perfect alignment of the controller coordinate system with the rotor coordinate system. Then, using :eq:`sm_model`, :eq:`flux_mapping_sm`, and :eq:`cc_flux`, the closed-loop system can be shown to be analogous to the induction machine case. This control design corresponds to the implementation in the :class:`motulator.drive.control.sm.CurrentController` class.
 
 .. rubric:: References
 
@@ -141,5 +141,3 @@ where :math:`\hat R_\mathrm{s} = 0` can be used in practice. Assume accurate par
 .. [#Awa2019] Awan, Saarakkala, Hinkkanen, "Flux-linkage-based current control of saturated synchronous motors," IEEE Trans. Ind. Appl. 2019, https://doi.org/10.1109/TIA.2019.2919258
 
 .. [#Hin2024] Hinkkanen,  Harnefors, Kukkola, "Fundamentals of Electric Machine Drives," lecture notes, 2024, https://doi.org/10.5281/zenodo.10609166
-
-
diff --git a/docs/source/control/drive/observers.rst b/docs/source/control/drive/observers.rst
index 59845c23..095084ab 100644
--- a/docs/source/control/drive/observers.rst
+++ b/docs/source/control/drive/observers.rst
@@ -15,10 +15,10 @@ The inverse-Γ model of an induction machine is considered (see :doc:`/model/dri
 
 .. math::
 	\frac{\mathrm{d} \boldsymbol{\psi}_\mathrm{R}}{\mathrm{d} t} + \mathrm{j}\omega_\mathrm{s}\boldsymbol{\psi}_\mathrm{R} &= \boldsymbol{u}_\mathrm{s} - R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t} - \mathrm{j} \omega_\mathrm{s}L_\sigma\boldsymbol{i}_\mathrm{s} \\
-    &= R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\alpha - \mathrm{j}\omega_\mathrm{m} \right)\boldsymbol{\psi}_\mathrm{R}
+    &= \underbrace{R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\alpha - \mathrm{j}\omega_\mathrm{m} \right)\boldsymbol{\psi}_\mathrm{R}}_{\boldsymbol{v}}
     :label: dpsiR
 
-where :math:`\omega_\mathrm{m}` is the electrical rotor angular speed.
+where :math:`\omega_\mathrm{m}` is the electrical rotor angular speed and :math:`\alpha = R_\mathrm{R}/L_\mathrm{M}` is the inverse rotor time constant.
 
 Observer Structure
 ^^^^^^^^^^^^^^^^^^
@@ -29,20 +29,20 @@ General Coordinates
 Based on :eq:`dpsiR`, a reduced-order observer can be formulated in a coordinate system rotating at :math:`\omega_\mathrm{s}`,
 
 .. math::
-    \frac{\mathrm{d} \hat{\boldsymbol{\psi}}_\mathrm{R}}{\mathrm{d} t} + \mathrm{j}\omega_\mathrm{s}\hat{\boldsymbol{\psi}}_\mathrm{R} = \boldsymbol{v} + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \boldsymbol{v}) + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \boldsymbol{v})^*
+    \frac{\mathrm{d} \hat{\boldsymbol{\psi}}_\mathrm{R}}{\mathrm{d} t} + \mathrm{j}\omega_\mathrm{s}\hat{\boldsymbol{\psi}}_\mathrm{R} = \hat{\boldsymbol{v}}_\mathrm{s} + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s}) + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s})^*
     :label: dhatpsiR
 
 where :math:`\boldsymbol{k}_1` and :math:`\boldsymbol{k}_2` are complex gains, the estimates are marked with the hat, and :math:`^*` marks the complex conjugate. The back-emf estimates are
 
 .. math::
-    \boldsymbol{v} &= \boldsymbol{u}_\mathrm{s} - R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t} - \mathrm{j} \omega_\mathrm{s}L_\sigma\boldsymbol{i}_\mathrm{s} \\
-	\hat{\boldsymbol{v}} &= R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\alpha - \mathrm{j}\hat{\omega}_\mathrm{m} \right)\hat{\boldsymbol{\psi}}_\mathrm{R}
+    \hat{\boldsymbol{v}}_\mathrm{s} &= \boldsymbol{u}_\mathrm{s} - \hat R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - \hat L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t} - \mathrm{j} \omega_\mathrm{s}\hat L_\sigma\boldsymbol{i}_\mathrm{s} \\
+	\hat{\boldsymbol{v}} &= \hat R_\mathrm{R}\boldsymbol{i}_\mathrm{s} - \left(\hat \alpha - \mathrm{j}\hat{\omega}_\mathrm{m} \right)\hat{\boldsymbol{\psi}}_\mathrm{R}
     :label: hatv
 
-In sensored drives, the rotor speed estimate is replaced with the measured speed,  :math:`\hat{\omega}_\mathrm{m} = \omega_\mathrm{m}`. If needed, the estimates for the stator flux linkage and the electromagnetic torque, respectively, are
+In sensored drives, the rotor speed estimate is replaced with the measured speed, :math:`\hat{\omega}_\mathrm{m} = \omega_\mathrm{m}`. If needed, the estimates for the stator flux linkage and the electromagnetic torque, respectively, are
 
 .. math::
-    \hat{\boldsymbol{\psi}}_\mathrm{s} &= \hat{\boldsymbol{\psi}}_\mathrm{R} + L_\sigma \boldsymbol{i}_\mathrm{s} \\
+    \hat{\boldsymbol{\psi}}_\mathrm{s} &= \hat{\boldsymbol{\psi}}_\mathrm{R} + \hat L_\sigma \boldsymbol{i}_\mathrm{s} \\
     \hat{\tau}_\mathrm{M} &= \frac{3n_\mathrm{p}}{2}\mathrm{Im}\left\{\boldsymbol{i}_\mathrm{s} \hat{\boldsymbol{\psi}}_\mathrm{R}^* \right\}
     :label: tauM
 
@@ -58,28 +58,28 @@ The flux estimate in stator coordinates can be expressed using the polar form, :
 
 .. math::
     \frac{\mathrm{d}\hat{\vartheta}_\mathrm{s}}{\mathrm{d} t} = \omega_\mathrm{s}
-    = \frac{\mathrm{Im} \{ \boldsymbol{v}' + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \boldsymbol{v}') + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \boldsymbol{v}')^* \} }{\hat{\psi}_\mathrm{R} + L_\sigma \mathrm{Re}\{(1 - \boldsymbol{k}_1)\boldsymbol{i}_\mathrm{s} + \boldsymbol{k}_2 \boldsymbol{i}_\mathrm{s}^* \}}
+    = \frac{\mathrm{Im} \{ \hat{\boldsymbol{v}}_\mathrm{s}' + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s}') + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s}')^* \} }{\hat{\psi}_\mathrm{R} + \hat L_\sigma \mathrm{Re}\{(1 - \boldsymbol{k}_1)\boldsymbol{i}_\mathrm{s} + \boldsymbol{k}_2 \boldsymbol{i}_\mathrm{s}^*\}}
     :label: hatws
 
 where
 
 .. math::
-    \boldsymbol{v}' = \boldsymbol{u}_\mathrm{s} - R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t}
+    \hat{\boldsymbol{v}}_\mathrm{s}' = \boldsymbol{u}_\mathrm{s} - \hat R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - \hat L_\sigma \frac{\mathrm{d} \boldsymbol{i}_\mathrm{s}}{\mathrm{d} t}
     :label: vp
 
 The flux magnitude dynamics are
 
 .. math::
     \frac{\mathrm{d} \hat{\psi}_\mathrm{R}}{\mathrm{d} t}
-    = \mathrm{Re}\{ \boldsymbol{v} + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \boldsymbol{v}) + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \boldsymbol{v})^* \}
+    = \mathrm{Re}\{ \hat{\boldsymbol{v}}_\mathrm{s} + \boldsymbol{k}_1(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s}) + \boldsymbol{k}_2(\hat{\boldsymbol{v}} - \hat{\boldsymbol{v}}_\mathrm{s})^* \}
     :label: dhatpsiR_abs
 
-Notice that the right-hand side of :eq:`hatws` is independent of :math:`\omega_\mathrm{s}`. Futhermore, in these coordinates, the condition :eq:`inherently` for an inherently sensorless observer reduces to :math:`\boldsymbol{k}_2 = \boldsymbol{k}_1`. This observer structure is implemented in the :class:`motulator.drive.control.im.Observer`, where simple forward-Euler discretization is used.
+Notice that the right-hand side of :eq:`hatws` is independent of :math:`\omega_\mathrm{s}`. Furthermore, in these coordinates, the condition :eq:`inherently` for an inherently sensorless observer reduces to :math:`\boldsymbol{k}_2 = \boldsymbol{k}_1`. This observer structure is implemented in the :class:`motulator.drive.control.im.Observer`, where simple forward-Euler discretization is used.
 
 Gain Selection
 ^^^^^^^^^^^^^^
 
-The estimation-error dynamics are obtained by subtracting :eq:`dhatpsiR` from :eq:`dpsiR`. The resulting system can be linearized for analysis and gain selection purposes. Using the rotor speed as an example, the small-signal deviation about the operating point is :math:`\Delta \omega_\mathrm{m} = \omega_\mathrm{m} - \omega_\mathrm{m0}`, where the subscript 0 refers to the operating point. Linearization of the estimation-error dynamics leads to [#Hin2010]_
+The estimation-error dynamics are obtained by subtracting :eq:`dhatpsiR` from :eq:`dpsiR`. The resulting system can be linearized for analysis and gain selection purposes. Accurate parameter estimates are assumed in the analysis. Using the rotor speed as an example, the small-signal deviation about the operating point is :math:`\Delta \omega_\mathrm{m} = \omega_\mathrm{m} - \omega_\mathrm{m0}`, where the subscript 0 refers to the operating point. Linearization of the estimation-error dynamics leads to [#Hin2010]_
 
 .. math::
 	\frac{\mathrm{d} \Delta\tilde{\boldsymbol{\psi}}_\mathrm{R}}{\mathrm{d} t} =  \boldsymbol{k}_1\Delta \tilde{\boldsymbol{v}} + \boldsymbol{k}_2\Delta \tilde{\boldsymbol{v}}^* - \mathrm{j}\omega_\mathrm{s0}\Delta\tilde{\boldsymbol{\psi}}_\mathrm{R}
@@ -103,7 +103,7 @@ Here, :math:`\boldsymbol{k}_2 = 0` and :math:`\Delta\tilde{\omega}_\mathrm{m} =
 It can be noticed that the closed-loop pole could be arbitrarily placed via the gain :math:`\boldsymbol{k}_1`. Well-damped estimation-error dynamics can be obtained, e.g., by choosing
 
 .. math::
-    \boldsymbol{k}_1 = 1 + \frac{g |\omega_\mathrm{m}|}{\alpha - \mathrm{j}\omega_\mathrm{m}}
+    \boldsymbol{k}_1 = 1 + \frac{g |\omega_\mathrm{m}|}{\hat \alpha - \mathrm{j}\omega_\mathrm{m}}
     :label: k1_sensored
 
 where :math:`g` is a unitless positive design parameter. The corresponding pole is located at :math:`s = -\alpha - g |\omega_\mathrm{m0}| - \mathrm{j}\omega_\mathrm{r0}`, where :math:`\omega_\mathrm{r0} = \omega_\mathrm{s0} - \omega_\mathrm{m0}` is the operating-point slip angular frequency.
@@ -131,10 +131,10 @@ yields an inherently sensorless observer, i.e., the rotor speed estimate :math:`
 where the gain components correspond to :math:`\boldsymbol{k}_1 = k_\mathrm{d} + \mathrm{j}k_\mathrm{q}`. It can be seen that the dynamics of the rotor speed are decoupled from the flux-estimation error dynamics. The decay rate :math:`\sigma` be assigned by choosing
 
 .. math::
-    \boldsymbol{k}_1 = \frac{\sigma}{\alpha - \mathrm{j}\hat\omega_\mathrm{m}}
+    \boldsymbol{k}_1 = \frac{\sigma}{\hat \alpha - \mathrm{j}\hat\omega_\mathrm{m}}
     :label: k1_sensorless
 
-which results in the characteristic polynomial :math:`D(s)=s^2 + 2\sigma s + \omega_\mathrm{s0}^2`. The decay rate can be selected as :math:`\sigma = \alpha/2 + \zeta_\infty|\hat{\omega}_\mathrm{m}|`, where :math:`\zeta_\infty` is the desired damping ratio at high speed. At zero stator frequency :math:`\omega_\mathrm{s0} = 0`, the poles are located at :math:`s = 0` and :math:`s = -\alpha`, which allows stable magnetizing and starting the machine.
+which results in the characteristic polynomial :math:`D(s)=s^2 + 2\sigma s + \omega_\mathrm{s0}^2`. The decay rate can be selected as :math:`\sigma = \hat \alpha/2 + \zeta_\infty|\hat{\omega}_\mathrm{m}|`, where :math:`\zeta_\infty` is the desired damping ratio at high speeds. At zero stator frequency :math:`\omega_\mathrm{s0} = 0`, the poles are located at :math:`s = 0` and :math:`s = -\alpha`, which allows stable magnetizing and starting the machine.
 
 .. rubric:: References
 
@@ -179,13 +179,13 @@ The observer is assumed to operate in estimated rotor coordinates, whose d-axis
 Since the stator current is measured, the observer is fundamentally corrected by means of the current estimation error. However, due to the saliency, it is more convenient to scale the current estimation error by the stator inductance, resulting in the flux linkage error
 
 .. math::
-	\boldsymbol{e} = \psi_\mathrm{f} + L_\mathrm{d} \mathrm{Re}\{ \boldsymbol{i}_\mathrm{s}\} + \mathrm{j} L_\mathrm{q} \mathrm{Im}\{\boldsymbol{i}_\mathrm{s} \} - \hat{\boldsymbol{\psi}}_\mathrm{s}
+	\boldsymbol{e} = \hat \psi_\mathrm{f} + \hat L_\mathrm{d} \mathrm{Re}\{ \boldsymbol{i}_\mathrm{s}\} + \mathrm{j} \hat L_\mathrm{q} \mathrm{Im}\{\boldsymbol{i}_\mathrm{s} \} - \hat{\boldsymbol{\psi}}_\mathrm{s}
     :label: e
 
 where :math:`\hat{\boldsymbol{\psi}}_\mathrm{s}` is the stator flux estimate. The flux linkage is estimated by
 
 .. math::
-    \frac{\mathrm{d} \hat{\boldsymbol{\psi}}_\mathrm{s}}{\mathrm{d} t} = \boldsymbol{u}_\mathrm{s} - R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - \mathrm{j}\omega_\mathrm{s}\hat{\boldsymbol{\psi}}_\mathrm{s} + \boldsymbol{k}_1 \boldsymbol{e} + \boldsymbol{k}_2 \boldsymbol{e}^*
+    \frac{\mathrm{d} \hat{\boldsymbol{\psi}}_\mathrm{s}}{\mathrm{d} t} = \boldsymbol{u}_\mathrm{s} - \hat R_\mathrm{s}\boldsymbol{i}_\mathrm{s} - \mathrm{j}\omega_\mathrm{s}\hat{\boldsymbol{\psi}}_\mathrm{s} + \boldsymbol{k}_1 \boldsymbol{e} + \boldsymbol{k}_2 \boldsymbol{e}^*
     :label: sm_flux_observer
 
 where :math:`\boldsymbol{k}_1` and :math:`\boldsymbol{k}_2` are gains (complex in a general case), the estimates are marked with the hat, and :math:`^*` marks the complex conjugate.
@@ -240,7 +240,7 @@ where :math:`\sigma` is the desired decay rate of the flux estimation error and
 allows to decouple the flux-estimation error dynamics from the rotor-position dynamics. By default, the decay rate is scheduled as
 
 .. math::
-    \sigma = \frac{R_\mathrm{s}}{4}\left(\frac{1}{L_\mathrm{d}} + \frac{1}{L_\mathrm{q}}\right) + \zeta_\infty |\hat{\omega}_\mathrm{m} |
+    \sigma = \frac{\hat R_\mathrm{s}}{4}\left(\frac{1}{\hat L_\mathrm{d}} + \frac{1}{\hat L_\mathrm{q}}\right) + \zeta_\infty |\hat{\omega}_\mathrm{m} |
     :label: sigma_sensorless
 
 where :math:`\zeta_\infty` is the desired damping ratio at high speed. At zero speed, :eq:`sigma_sensorless` places one pole at :math:`s = 0` and another at :math:`s = -(R_\mathrm{s}/2)(1/L_\mathrm{d} + 1/L_\mathrm{q})`.
diff --git a/docs/source/control/grid/current_ctrl.rst b/docs/source/control/grid/current_ctrl.rst
index 8d279e8e..9306a3b7 100644
--- a/docs/source/control/grid/current_ctrl.rst
+++ b/docs/source/control/grid/current_ctrl.rst
@@ -1,7 +1,7 @@
 Current Control
 ===============
 
-Synchronous-frame two-degrees-of-freedom (2DOF) proportional-integral (PI) current control can be used in grid converters [#Har2015]_. This control structure allows to compensate for the cross-coupling originating from rotating coordinates as well as to improve disturbance rejection. A 2DOF PI current controller is available in the :class:`motulator.grid.control.CurrentController` class, whose base class is :class:`motulator.common.control.ComplexPIController`. 
+Synchronous-frame two-degrees-of-freedom (2DOF) proportional-integral (PI) current control can be used in grid converters [#Har2015]_. This control structure allows to compensate for the cross-coupling originating from rotating coordinates as well as to improve disturbance rejection. A 2DOF PI current controller is available in the :class:`motulator.grid.control.CurrentController` class, whose base class is :class:`motulator.common.control.ComplexPIController`.
 
 .. note:: This controller design assumes an L filter, but it can also be applied with LCL filters (see the :doc:`/grid_examples/grid_following/plot_gfl_lcl_10kva` example). If LCL-resonance damping and very low sampling frequencies are needed, the controller could be designed directly in the discrete-time domain taking the LCL filter dynamics into account [#Rah2021]_.
 
@@ -12,13 +12,13 @@ The design of synchronous-frame 2DOF PI current control is considered in the con
 
 .. math::
 	\frac{\mathrm{d} \boldsymbol{u}_\mathrm{i}}{\mathrm{d} t} &= (\boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{c}\boldsymbol{k}_\mathrm{t} )\left(\boldsymbol{i}_\mathrm{c,ref} - \boldsymbol{i}_\mathrm{c}\right) \\
-    \boldsymbol{u}_\mathrm{c,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{i}_\mathrm{c,ref} - \boldsymbol{k}_\mathrm{p}\boldsymbol{i}_\mathrm{c} + \boldsymbol{u}_\mathrm{i} 
+    \boldsymbol{u}_\mathrm{c,ref} &= \boldsymbol{k}_\mathrm{t}\boldsymbol{i}_\mathrm{c,ref} - \boldsymbol{k}_\mathrm{p}\boldsymbol{i}_\mathrm{c} + \boldsymbol{u}_\mathrm{i}
     :label: grid_cc
 
-where :math:`\boldsymbol{u}_\mathrm{c,ref}` is the output of the controller, i.e., the converter voltage reference, :math:`\boldsymbol{i}_\mathrm{c}` is the measured converter current, :math:`\boldsymbol{i}_\mathrm{c,ref}` is the converter current reference, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\omega_\mathrm{c}` is the angular speed of the coordinate system. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain. 
+where :math:`\boldsymbol{u}_\mathrm{c,ref}` is the output of the controller, i.e., the converter voltage reference, :math:`\boldsymbol{i}_\mathrm{c}` is the measured converter current, :math:`\boldsymbol{i}_\mathrm{c,ref}` is the converter current reference, :math:`\boldsymbol{u}_\mathrm{i}` is the the integral state, and :math:`\omega_\mathrm{c}` is the angular speed of the coordinate system. Furthermore, :math:`\boldsymbol{k}_\mathrm{t}` is the reference-feedforward gain, :math:`\boldsymbol{k}_\mathrm{p}` is the proportional gain, and :math:`\boldsymbol{k}_\mathrm{i}` is the integral gain.
 
-Closed-Loop System 
-^^^^^^^^^^^^^^^^^^
+Closed-Loop System
+------------------
 
 Consider the grid model in synchronous coordinates
 
@@ -34,23 +34,23 @@ where :math:`\boldsymbol{u}_\mathrm{c}` is the converter output voltage, :math:`
 The disturbance rejection depends on the closed-loop admittance
 
 .. math::
-    \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L s^2 + (\boldsymbol{k}_\mathrm{p} + \mathrm{j}\omega_\mathrm{c} L) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{c} \boldsymbol{k}_\mathrm{t}} 
+    \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L s^2 + (\boldsymbol{k}_\mathrm{p} + \mathrm{j}\omega_\mathrm{c} L) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{c} \boldsymbol{k}_\mathrm{t}}
     :label: Yc_grid
 
 The closed-loop poles can be arbitrarily placed by means of the gains. The reference-tracking transfer function is
 
 .. math::
-	\boldsymbol{G}_\mathrm{c}(s) = \frac{(s + \mathrm{j}\omega_\mathrm{c}) \boldsymbol{k}_\mathrm{t} + \boldsymbol{k}_\mathrm{i} }{L s^2 + (\boldsymbol{k}_\mathrm{p} + \mathrm{j}\omega_\mathrm{c} L) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{c} \boldsymbol{k}_\mathrm{t}}     
+	\boldsymbol{G}_\mathrm{c}(s) = \frac{(s + \mathrm{j}\omega_\mathrm{c}) \boldsymbol{k}_\mathrm{t} + \boldsymbol{k}_\mathrm{i} }{L s^2 + (\boldsymbol{k}_\mathrm{p} + \mathrm{j}\omega_\mathrm{c} L) s + \boldsymbol{k}_\mathrm{i} + \mathrm{j}\omega_\mathrm{c} \boldsymbol{k}_\mathrm{t}}
     :label: Gc_grid
 
 whose zero can be placed by means of the reference-feedforward gain :math:`\boldsymbol{k}_\mathrm{t}`.
 
 Gain Selection
-^^^^^^^^^^^^^^
+--------------
 
 Consider the gains
 
-.. math::                
+.. math::
     \boldsymbol{k}_\mathrm{p} = 2\alpha_\mathrm{c} \hat L \qquad\qquad
     \boldsymbol{k}_\mathrm{i} = \alpha_\mathrm{c}^2\hat L  \qquad \qquad
     \boldsymbol{k}_\mathrm{t} = \alpha_\mathrm{c} \hat L
@@ -60,9 +60,9 @@ where :math:`\hat L` is the inductance estimate. Assuming accurate parameter est
 .. math::
     \boldsymbol{Y}_\mathrm{c}(s) = \frac{s}{L (s + \alpha_\mathrm{c})(s + \alpha_\mathrm{c} + \mathrm{j}\omega_\mathrm{c})}
     \qquad\qquad
-    \boldsymbol{G}_\mathrm{c}(s) = \frac{\alpha_\mathrm{c}}{s + \alpha_\mathrm{c}} 
+    \boldsymbol{G}_\mathrm{c}(s) = \frac{\alpha_\mathrm{c}}{s + \alpha_\mathrm{c}}
 
-It can be seen that this design results in the first-order reference-tracking dynamics. Furthermore, one pole is placed at the real axis at :math:`s=-\alpha_\mathrm{c}` and another pole at :math:`s= -\alpha_\mathrm{c} - \mathrm{j}\omega_\mathrm{c}`. This gain selection is used in the :class:`motulator.grid.control.CurrentController` class. 
+It can be seen that this design results in the first-order reference-tracking dynamics. Furthermore, one pole is placed at the real axis at :math:`s=-\alpha_\mathrm{c}` and another pole at :math:`s= -\alpha_\mathrm{c} - \mathrm{j}\omega_\mathrm{c}`. This gain selection is used in the :class:`motulator.grid.control.CurrentController` class.
 
 The converter output voltage is limited in practice due to the limited DC-bus voltage of the converter. Consequently, the realized (limited) voltage reference is
 
@@ -80,6 +80,3 @@ where :math:`\mathrm{sat}(\cdot)` is the saturation function. The limited voltag
 .. [#Har1998] Harnefors, Nee, "Model-based current control of AC machines using the internal model control method," IEEE Trans. Ind. Appl., 1998, https://doi.org/10.1109/28.658735
 
 .. [#Bri1999] Briz del Blanco, Degner, Lorenz, “Dynamic analysis of current regulators for AC motors using complex vectors,” IEEE Trans. Ind. Appl., 1999, https://doi.org/10.1109/28.806058
-
-
-
diff --git a/docs/source/control/grid/dc_voltage_ctrl.rst b/docs/source/control/grid/dc_voltage_ctrl.rst
index 23bc2ed6..d9cbd2cf 100644
--- a/docs/source/control/grid/dc_voltage_ctrl.rst
+++ b/docs/source/control/grid/dc_voltage_ctrl.rst
@@ -1,7 +1,7 @@
 DC-Bus Voltage Control
 ======================
 
-A PI controller for the DC-bus voltage is implemented in the :class:`motulator.grid.control.DCBusVoltageController` class, whose base class is :class:`motulator.common.control.PIController`. In the following, the tuning of the DC-bus voltage controller is considered in the continuous-time domain. 
+A PI controller for the DC-bus voltage is implemented in the :class:`motulator.grid.control.DCBusVoltageController` class, whose base class is :class:`motulator.common.control.PIController`. In the following, the tuning of the DC-bus voltage controller is considered in the continuous-time domain.
 
 Controller Structure
 --------------------
@@ -10,16 +10,16 @@ In order to have a linear closed-loop system, it is convenient to use the DC-bus
 
 .. math::
 	\hat W_\mathrm{dc} = \frac{1}{2}\hat{C}_\mathrm{dc} u_\mathrm{dc}^2 \qquad
-    W_\mathrm{dc,ref} = \frac{1}{2}\hat{C}_\mathrm{dc} u_\mathrm{dc,ref}^2 
-    :label: dc_cap_energy 
+    W_\mathrm{dc,ref} = \frac{1}{2}\hat{C}_\mathrm{dc} u_\mathrm{dc,ref}^2
+    :label: dc_cap_energy
 
-where :math:`\hat{C}_\mathrm{dc}` is the DC-bus capacitance estimate, :math:`u_\mathrm{dc}` is the measured DC-bus voltage, and :math:`u_\mathrm{dc,ref}` is the DC-bus voltage reference. Using these variables, the output power reference :math:`p_\mathrm{c,ref}` for the converter is obtained from a PI controller, 
+where :math:`\hat{C}_\mathrm{dc}` is the DC-bus capacitance estimate, :math:`u_\mathrm{dc}` is the measured DC-bus voltage, and :math:`u_\mathrm{dc,ref}` is the DC-bus voltage reference. Using these variables, the output power reference :math:`p_\mathrm{c,ref}` for the converter is obtained from a PI controller,
 
 .. math::
     p_\mathrm{c,ref} = -k_\mathrm{p} (W_\mathrm{dc,ref} - \hat{W}_\mathrm{dc}) - \int k_\mathrm{i} (W_\mathrm{dc,ref} - \hat{W}_\mathrm{dc}) \mathrm{d} t
     :label: dc_voltage_ctrl
 
-where :math:`k_\mathrm{p}` is the proportional gain and :math:`k_\mathrm{i}` is the integral gain. The negative signs appear since the positive converter output power decreases the DC-bus capacitor energy, see :eq:`capacitor`. 
+where :math:`k_\mathrm{p}` is the proportional gain and :math:`k_\mathrm{i}` is the integral gain. The negative signs appear since the positive converter output power decreases the DC-bus capacitor energy, see :eq:`capacitor`.
 
 Closed-Loop System
 ------------------
@@ -38,7 +38,7 @@ In the Laplace domain, the closed-loop system resulting from :eq:`dc_voltage_ctr
     W_\mathrm{dc}(s) = \frac{k_\mathrm{p} s + k_\mathrm{i}}{s^2 + k_\mathrm{p} s + k_\mathrm{i}} W_\mathrm{dc,ref}(s) + \frac{s}{s^2 + k_\mathrm{p} s + k_\mathrm{i}} p_\mathrm{dc}(s)
     :label: dc_voltage_ctrl_closed_loop
 
-where it can be seen that :math:`W_\mathrm{dc} = W_\mathrm{dc,ref}` holds in the steady state. Furthermore, it can be shown that also :math:`u_\mathrm{dc} = u_\mathrm{dc,ref}` holds in the steady state (independently of the errors in the capacitance estimate :math:`\hat{C}`, since both reference and measured values in :eq:`dc_cap_energy` are scaled by the same estimate). 
+where it can be seen that :math:`W_\mathrm{dc} = W_\mathrm{dc,ref}` holds in the steady state. Furthermore, it can be shown that also :math:`u_\mathrm{dc} = u_\mathrm{dc,ref}` holds in the steady state (independently of the errors in the capacitance estimate :math:`\hat{C}`, since both reference and measured values in :eq:`dc_cap_energy` are scaled by the same estimate).
 
 Gain Selection
 --------------
@@ -47,9 +47,9 @@ Based on :eq:`dc_voltage_ctrl_closed_loop`, the gain selection
 
 .. math::
     k_\mathrm{p} = 2\alpha_\mathrm{dc} \qquad
-    k_\mathrm{i} = \alpha_\mathrm{dc}^2  
+    k_\mathrm{i} = \alpha_\mathrm{dc}^2
 
-results in the double real pole at :math:`s = -\alpha_\mathrm{dc}`. The closed-loop bandwidth is approximately :math:`\alpha_\mathrm{dc}`. 
+results in the double real pole at :math:`s = -\alpha_\mathrm{dc}`. The closed-loop bandwidth is approximately :math:`\alpha_\mathrm{dc}`.
 
 .. rubric:: References
 
diff --git a/docs/source/control/grid/grid-forming.rst b/docs/source/control/grid/grid-forming.rst
new file mode 100644
index 00000000..939e90c6
--- /dev/null
+++ b/docs/source/control/grid/grid-forming.rst
@@ -0,0 +1,92 @@
+Disturbance-Observer-Based Grid-Forming Control
+===============================================
+
+In these notes, disturbance-observer-based grid-forming control is discussed [#Nur2024]_. This control method is available in the :class:`motulator.grid.control.ObserverBasedGridFormingControl` class. As can be realized by comparing the examples :doc:`/grid_examples/grid_forming/plot_gfm_rfpsc_13kva` and :doc:`/grid_examples/grid_forming/plot_gfm_obs_13kva`, the observer-based grid-forming control can be configured to behave practically identically with reference-feedforward power-synchronization control (PSC) [#Har2020]_. Compared to reference-feedforward PSC, disturbance-observer-based grid-forming control is easier to analyze and extend with different operating modes (including grid-following and transparent current-control modes). The order of these two control methods is the same and their computational complexity is comparable.
+
+System Model
+------------
+
+First, the system is modeled in general coordinates, whose angle with respect to the stationary coordinates is :math:`\vartheta_\mathrm{c}` and the angular speed is :math:`\omega_\mathrm{c} = \mathrm{d}\vartheta_\mathrm{c}/\mathrm{d} t`. The dynamics of the inductor current :math:`\boldsymbol{i}_\mathrm{c}` and the grid voltage :math:`\boldsymbol{u}_\mathrm{g}` are modeled as
+
+.. math::
+    L\frac{\mathrm{d} \boldsymbol{i}_\mathrm{c}}{\mathrm{d} t} &= \boldsymbol{u}_\mathrm{c} - \boldsymbol{u}_\mathrm{g} - \mathrm{j}\omega_\mathrm{c} L \boldsymbol{i}_\mathrm{c} \\
+    \frac{\mathrm{d}\boldsymbol{u}_\mathrm{g}}{\mathrm{d} t} &= \mathrm{j}(\omega_\mathrm{g} - \omega_\mathrm{c})\boldsymbol{u}_\mathrm{g}
+    :label: system_model1
+
+where :math:`L` is the inductance and :math:`\omega_\mathrm{g}` is the grid angular frequency. The grid voltage :math:`\boldsymbol{u}_\mathrm{g}` is modeled as a disturbance; this same disturbance model for the grid voltage is also used in the development of :doc:`/control/grid/pll`. The active power fed to the grid is nonlinear in the state variables
+
+.. math::
+    p_\mathrm{g} = \frac{3}{2}\mathrm{Re}\{\boldsymbol{u}_\mathrm{g}\boldsymbol{i}_\mathrm{c}^*\}
+    :label: power1
+
+For the purpose of grid-forming control, we may define the quasi-static converter voltage as an output variable as a function of the state variables,
+
+.. math::
+    \boldsymbol{v}_\mathrm{c} = \boldsymbol{u}_\mathrm{g} + \mathrm{j}\omega_\mathrm{g} L \boldsymbol{i}_\mathrm{c}
+    :label: quasi_static_converter_voltage
+
+This voltage is controllable (unlike the grid voltage which is the disturbance) and it equals the converter output voltage :math:`\boldsymbol{u}_\mathrm{c}` in the steady state. The disturbance-observer-based PSC regulates the voltage magnitude :math:`|\boldsymbol{v}_\mathrm{c}|` and the power :math:`p_\mathrm{g}` fed to grid.
+
+Control System
+--------------
+
+The control system consists of a disturbance observer and a control law. The disturbance observer estimates the quasi-static converter voltage and the power fed to the grid. The control law generates the converter voltage reference based on the estimated grid voltage and power. The disturbance observer also provides the integral action for the control law. This control structure allows seamless switching between grid-forming and grid-following modes [#Nur2024]_. For simplicity, only the grid-forming mode is considered in the following.
+
+Disturbance Observer
+^^^^^^^^^^^^^^^^^^^^
+
+Based on :eq:`system_model1`--:eq:`quasi_static_converter_voltage`, a disturbance observer for the grid voltage can be formed [#Nur2024]_, [#Kuk2021]_, [#Fra1997]_
+
+.. math::
+    \frac{\mathrm{d} \hat{\boldsymbol{u}}_\mathrm{g}}{\mathrm{d} t} &= \mathrm{j} (\hat{\omega}_\mathrm{g} - \omega_\mathrm{c})\hat{\boldsymbol{u}}_\mathrm{g} + \alpha_\mathrm{o}\left(\boldsymbol{u}_\mathrm{c} - \hat L \frac{\mathrm{d} \boldsymbol{i}_\mathrm{c}}{\mathrm{d} t} - \mathrm{j} \omega_\mathrm{c} \hat L \boldsymbol{i}_\mathrm{c} - \hat{\boldsymbol{u}}_\mathrm{g} \right) \\
+    \hat{\boldsymbol{v}}_\mathrm{c} &= \hat{\boldsymbol{u}}_\mathrm{g} + \mathrm{j}\hat{\omega}_\mathrm{g} \hat L \boldsymbol{i}_\mathrm{c} \\
+    \hat p_\mathrm{g} &= \frac{3}{2}\mathrm{Re}\{\hat{\boldsymbol{v}}_\mathrm{c}\boldsymbol{i}_\mathrm{c}^*\}
+    :label: disturbance_observer_gfm
+
+where :math:`\hat \omega_\mathrm{g}` is the nominal grid angular frequency and estimates are marked with hat.
+
+.. note:: A conventional phase-locked loop (PLL) can be expressed in the same disturbance observer framework, see the :doc:`/control/grid/pll` notes. It can be realized that the measured grid voltage :math:`\boldsymbol{u}_\mathrm{g}` used in the conventional PLL is replaced by its converter-voltage-based estimate :math:`\boldsymbol{u}_\mathrm{c} - \hat L (\mathrm{d} \boldsymbol{i}_\mathrm{c}/\mathrm{d} t) - \mathrm{j} \omega_\mathrm{c} \hat L \boldsymbol{i}_\mathrm{c}` in the disturbance observer :eq:`disturbance_observer_gfm`.
+
+Control Law
+^^^^^^^^^^^
+
+A nonlinear state feedback law is used
+
+.. math::
+    \boldsymbol{u}_\mathrm{c,ref} = \hat{\boldsymbol{v}}_\mathrm{c} + \boldsymbol{k}_\mathrm{p} (p_\mathrm{g,ref} - \hat p_\mathrm{g}) + \boldsymbol{k}_\mathrm{v} (v_\mathrm{c,ref} - |\hat{\boldsymbol{v}}_\mathrm{c}|)
+    :label: control_law_gfm
+
+where :math:`p_\mathrm{g,ref}` is the active power reference, :math:`v_\mathrm{c,ref}` is the converter voltage magnitude reference, and :math:`\boldsymbol{k}_\mathrm{p}` and :math:`\boldsymbol{k}_\mathrm{v}` are the complex gains for the active-power and converter-voltage-magnitude channels, respectively. To ensure robust operation, similar to the reference-feedforward PSC, the complex gains can be selected as
+
+.. math::
+    \boldsymbol{k}_\mathrm{p} = \frac{R_\mathrm{a}}{v_\mathrm{c,ref}} \frac{\hat{\boldsymbol{v}}_\mathrm{c}}{|\hat{\boldsymbol{v}}_\mathrm{c}|} \qquad
+    \boldsymbol{k}_\mathrm{v} = (1 - \mathrm{j} k_\mathrm{v}) \frac{\hat{\boldsymbol{v}}_\mathrm{c}}{|\hat{\boldsymbol{v}}_\mathrm{c}|}
+    :label: gain_selection_gfm
+
+where the gains :math:`R_\mathrm{a} = 0.2` p.u. and :math:`k_\mathrm{v} = \alpha_\mathrm{o}/\omega_\mathrm{g}` can be used.
+
+Implementation Aspects
+^^^^^^^^^^^^^^^^^^^^^^
+
+To avoid the derivate on the right-hand side of :eq:`disturbance_observer_gfm`, a new state variable :math:`\hat{\boldsymbol{u}}_\mathrm{g}' = \hat{\boldsymbol{u}}_\mathrm{g} + \alpha_\mathrm{o} \hat L \boldsymbol{i}_\mathrm{c}` can be introduced [#Fra1997]_. Furthermore, the coordinate system for the implementation can be chosen freely. The simplest choice is to use the nominal grid frequency as the coordinate system frequency, :math:`\omega_\mathrm{c} = \hat \omega_\mathrm{g}`. Using these design choices, the whole control system consisting of the disturbance observer :eq:`disturbance_observer_gfm` and the control law :eq:`control_law_gfm` in the state-space form reduces to [#Nur2024]_
+
+.. math::
+    \frac{\mathrm{d} \hat{\boldsymbol{u}}'_\mathrm{g}}{\mathrm{d} t} &= \alpha_\mathrm{o} (\boldsymbol{u}_\mathrm{c,ref} - \hat{\boldsymbol{v}}_\mathrm{c} ) \\
+    \hat{\boldsymbol{v}}_\mathrm{c} &= \hat{\boldsymbol{u}}_\mathrm{g}' - (\alpha_\mathrm{o} - \mathrm{j}\hat{\omega}_\mathrm{g}) \hat L \boldsymbol{i}_\mathrm{c} \\
+    \hat p_\mathrm{g} &= \frac{3}{2}\mathrm{Re}\{\hat{\boldsymbol{v}}_\mathrm{c}\boldsymbol{i}_\mathrm{c}^*\} \\
+    \boldsymbol{u}_\mathrm{c,ref} &= \hat{\boldsymbol{v}}_\mathrm{c} + \boldsymbol{k}_\mathrm{p} (p_\mathrm{g,ref} - \hat p_\mathrm{g}) + \boldsymbol{k}_\mathrm{v} (v_\mathrm{c,ref} - |\hat{\boldsymbol{v}}_\mathrm{c}|)
+    :label: control_system_gfm
+
+where the gains can be selected according to :eq:`gain_selection_gfm` and the converter voltage appearing in the observer has been replaced with its reference. Various control modes could be easily incorporated into the control system :eq:`control_system_gfm`, simply by changing the feedback correction terms of the control law [#Nur2024]_. The switching between the modes is seamless since the control law does not have memory, but the integral action is provided by the disturbance observer (in addition to synchronization).
+
+The control system implemented in the :class:`motulator.grid.control.ObserverBasedGridFormingControl` class corresponds to :eq:`control_system_gfm`. In the example implementation, a transparent current-control mode is implemented. In the grid-forming mode, the observer bandwidth :math:`\alpha_\mathrm{o} = 1` p.u. can be used. Furthermore, the inductance estimate can be set close to the lowest expected inductance value, e.g., :math:`\hat L = 0.15` p.u. Using this configuration, the robust performance from strong grids to very weak grids can be achieved. This grid-forming control method can also be used with LCL filters, similarly to reference-feedforward PSC.
+
+.. rubric:: References
+
+.. [#Nur2024] Nurminen, Mourouvin, Hinkkanen, Kukkola, "Multifunctional grid-forming converter control based on a disturbance observer, "IEEE Trans. Power Electron., 2024, https://doi.org/10.1109/TPEL.2024.3433503
+
+.. [#Har2020] Harnefors, Rahman, Hinkkanen, Routimo, "Reference-feedforward power-synchronization control," IEEE Trans. Power Electron., 2020, https://doi.org/10.1109/TPEL.2020.2970991
+
+.. [#Kuk2021] Kukkola, Routimo, Hinkkanen, Harnefors, "A voltage-sensorless controller for grid converters," IEEE PES ISGT Europe, 2021, https://doi.org/10.1109/ISGTEurope52324.2021.9640206
+
+.. [#Fra1997] Franklin, Powell, Workman, "Digital Control of Dynamic Systems," 3rd ed., Menlo Park, CA: Addison-Wesley, 1997
diff --git a/docs/source/control/grid/index.rst b/docs/source/control/grid/index.rst
index c91e97d3..4e4a3f0e 100644
--- a/docs/source/control/grid/index.rst
+++ b/docs/source/control/grid/index.rst
@@ -2,7 +2,7 @@
 Grid Converters
 ***************
 
-Design notes for selected control methods for grid-connected converters are provided in this section. The aim of these notes is to link the implemented methods to the theory and to provide a reference for the implementation. Further details are available in the references provided. 
+Design notes for selected control methods for grid-connected converters are provided in this section. The aim of these notes is to link the implemented methods to the theory and to provide a reference for the implementation. Further details are available in the references provided.
 
 .. toctree::
    :titlesonly:
@@ -10,4 +10,5 @@ Design notes for selected control methods for grid-connected converters are prov
 
    dc_voltage_ctrl
    current_ctrl
-   synchronization
+   Phase-Locked Loop <pll>
+   Grid-Forming Control <grid-forming>
diff --git a/docs/source/control/grid/synchronization.rst b/docs/source/control/grid/pll.rst
similarity index 78%
rename from docs/source/control/grid/synchronization.rst
rename to docs/source/control/grid/pll.rst
index dd297e0d..31986724 100644
--- a/docs/source/control/grid/synchronization.rst
+++ b/docs/source/control/grid/pll.rst
@@ -1,7 +1,7 @@
-Disturbance-Observer-Based PLL
-==============================
+Phase-Locked Loop
+=================
 
-A phase-locked loop (PLL) is commonly used in grid-following converters to synchronize the converter output with the grid [#Kau1997]_. Here, we represent the PLL using the disturbance observer structure [#Fra1997]_, which may be simpler to extend than the classical PLL. Furthermore, this structure allows to see links to synchronization methods used in grid-forming converters, see [#Nur2024]_. 
+A phase-locked loop (PLL) is commonly used in grid-following converters to synchronize the converter output with the grid [#Kau1997]_. Here, we represent the PLL using the disturbance observer structure [#Fra1997]_, which may be simpler to extend than the classical PLL. Furthermore, this structure allows to see links to synchronization methods used in grid-forming converters, see [#Nur2024]_.
 
 Disturbance Model
 -----------------
@@ -9,24 +9,33 @@ Disturbance Model
 Consider the positive-sequence voltage at the point of common coupling (PCC) in general coordinates, rotating at the angular speed :math:`\omega_\mathrm{c}`. The dynamics of the PCC voltage :math:`\boldsymbol{u}_\mathrm{g}` can be modeled using a disturbance model as
 
 .. math::
-    \frac{\mathrm{d}\boldsymbol{u}_\mathrm{g}}{\mathrm{d}t} = \mathrm{j}(\omega_\mathrm{g} - \omega_\mathrm{c})\boldsymbol{u}_\mathrm{g} 
+    \frac{\mathrm{d}\boldsymbol{u}_\mathrm{g}}{\mathrm{d}t} = \mathrm{j}(\omega_\mathrm{g} - \omega_\mathrm{c})\boldsymbol{u}_\mathrm{g}
     :label: disturbance
 
-where :math:`\omega_\mathrm{g}` is the grid angular frequency. If needed, this disturbance model could be extended, e.g., with a negative-sequence component, allowing to design the PLL for unbalanced grids. 
+where :math:`\omega_\mathrm{g}` is the grid angular frequency. If needed, this disturbance model could be extended, e.g., with a negative-sequence component, allowing to design the PLL for unbalanced grids.
+
+.. note:: The disturbance model in :eq:`disturbance` can be equivalently expressed in polar form as
+
+    .. math::
+        \frac{\mathrm{d} u_\mathrm{g}}{\mathrm{d} t} &= 0 \\
+        \frac{\mathrm{d} \delta_\mathrm{g}}{\mathrm{d} t} &= \omega_\mathrm{g} - \omega_\mathrm{c} \\
+        \boldsymbol{u}_\mathrm{g} &= u_\mathrm{g} \mathrm{e}^{\mathrm{j}\delta_\mathrm{g}}
+
+    where :math:`\delta_\mathrm{g} = \vartheta_\mathrm{g} - \vartheta_\mathrm{c}` is the angle of the grid voltage in general coordinates and the last equation is the output equation.
 
 PLL in General Coordinates
 --------------------------
 
-Based on :eq:`disturbance`, the disturbance observer containing the regular PLL functionality can be formulated as  
+Based on :eq:`disturbance`, the disturbance observer containing the regular PLL functionality can be formulated as
 
 .. math::
-    \frac{\mathrm{d}\hat{\boldsymbol{u}}_\mathrm{g}}{\mathrm{d}t} = \mathrm{j}(\hat{\omega}_\mathrm{g} - \omega_\mathrm{c})\hat{\boldsymbol{u}}_\mathrm{g} + \alpha_\mathrm{g} (\boldsymbol{u}_\mathrm{g} - \hat{\boldsymbol{u}}_\mathrm{g} ) 
+    \frac{\mathrm{d}\hat{\boldsymbol{u}}_\mathrm{g}}{\mathrm{d}t} = \mathrm{j}(\hat{\omega}_\mathrm{g} - \omega_\mathrm{c})\hat{\boldsymbol{u}}_\mathrm{g} + \alpha_\mathrm{o} (\boldsymbol{u}_\mathrm{g} - \hat{\boldsymbol{u}}_\mathrm{g} )
     :label: pll
 
-where :math:`\hat{\boldsymbol{u}}_\mathrm{g}` is the estimated PCC voltage, :math:`\hat{\omega}_\mathrm{g}` is the grid angular frequency estimate (either constant corresponding to the nominal value or dynamic from grid-frequency tracking), and :math:`\alpha_\mathrm{g}` is the bandwidth. If needed, the disturbance observer can be extended with the frequency tracking as
+where :math:`\hat{\boldsymbol{u}}_\mathrm{g}` is the estimated PCC voltage, :math:`\hat{\omega}_\mathrm{g}` is the grid angular frequency estimate (either constant corresponding to the nominal value or dynamic from grid-frequency tracking), and :math:`\alpha_\mathrm{o}` is the bandwidth. If needed, the disturbance observer can be extended with the frequency tracking as
 
 .. math::
-    \frac{\mathrm{d}\hat{\omega}_\mathrm{g}}{\mathrm{d}t} = k_\omega\mathrm{Im}\left\{ \frac{\boldsymbol{u}_\mathrm{g} - \hat{\boldsymbol{u}}_\mathrm{g}}{\hat{\boldsymbol{u}}_\mathrm{g}} \right\} 
+    \frac{\mathrm{d}\hat{\omega}_\mathrm{g}}{\mathrm{d}t} = k_\omega\mathrm{Im}\left\{ \frac{\boldsymbol{u}_\mathrm{g} - \hat{\boldsymbol{u}}_\mathrm{g}}{\hat{\boldsymbol{u}}_\mathrm{g}} \right\}
     :label: frequency_tracking
 
 where :math:`k_\omega` is the frequency-tracking gain. Notice that :eq:`pll` and :eq:`frequency_tracking` are both driven by the estimation error :math:`\boldsymbol{u}_\mathrm{g} - \hat{\boldsymbol{u}}_\mathrm{g}` of the PCC voltage.
@@ -37,8 +46,8 @@ PLL in Estimated PCC Voltage Coordinates
 The disturbance observer :eq:`pll` can be equivalently expressed in estimated PCC voltage coordinates, where :math:`\hat{\boldsymbol{u}}_\mathrm{g} = \hat{u}_\mathrm{g}` and :math:`\omega_\mathrm{c} = \mathrm{d} \hat{\vartheta}_\mathrm{g}/ \mathrm{d} t`, resulting in
 
 .. math::
-    \frac{\mathrm{d}\hat{u}_\mathrm{g}}{\mathrm{d}t} &= \alpha_\mathrm{g} \left(\mathrm{Re}\{\boldsymbol{u}_\mathrm{g}\} - \hat{u}_\mathrm{g} \right) \\
-    \frac{\mathrm{d} \hat{\vartheta}_\mathrm{g}}{\mathrm{d} t} &= \hat{\omega}_\mathrm{g} + \frac{\alpha_\mathrm{g}}{\hat{u}_\mathrm{g}}\mathrm{Im}\{ \boldsymbol{u}_\mathrm{g} \} = \omega_\mathrm{c}
+    \frac{\mathrm{d}\hat{u}_\mathrm{g}}{\mathrm{d}t} &= \alpha_\mathrm{o} \left(\mathrm{Re}\{\boldsymbol{u}_\mathrm{g}\} - \hat{u}_\mathrm{g} \right) \\
+    \frac{\mathrm{d} \hat{\vartheta}_\mathrm{g}}{\mathrm{d} t} &= \hat{\omega}_\mathrm{g} + \frac{\alpha_\mathrm{o}}{\hat{u}_\mathrm{g}}\mathrm{Im}\{ \boldsymbol{u}_\mathrm{g} \} = \omega_\mathrm{c}
     :label: pll_polar
 
 It can be seen that the first equation in :eq:`pll_polar` is low-pass filtering of the measured PCC voltage magnitude and the second equation is the conventional angle-tracking PLL. In these coordinates, the frequency tracking :eq:`frequency_tracking` reduces to
@@ -47,7 +56,7 @@ It can be seen that the first equation in :eq:`pll_polar` is low-pass filtering
     \frac{\mathrm{d} \hat{\omega}_\mathrm{g}}{\mathrm{d} t} = \frac{k_\omega}{\hat{u}_\mathrm{g} } \mathrm{Im}\{ \boldsymbol{u}_\mathrm{g} \}
     :label: frequency_tracking_polar
 
-It can be noticed that the disturbance observer with the frequency tracking equals the conventional frequency-adaptive PLL [#Kau1997]_, with the additional feature of low-pass filtering the PCC voltage magnitude. The low-pass filtered PCC voltage can be used as a feedforward term in current control [#Har2021]_. 
+It can be noticed that the disturbance observer with the frequency tracking equals the conventional frequency-adaptive PLL [#Kau1997]_, with the additional feature of low-pass filtering the PCC voltage magnitude. The low-pass filtered PCC voltage can be used as a feedforward term in current control [#Har2021]_.
 
 Linearized Closed-Loop System
 -----------------------------
@@ -55,37 +64,37 @@ Linearized Closed-Loop System
 The estimation-error dynamics are analyzed by means of linearization. Using the PCC voltage as an example, the small-signal deviation about the operating point is denoted by :math:`\Delta \boldsymbol{u}_\mathrm{g} = \boldsymbol{u}_\mathrm{g} - \boldsymbol{u}_\mathrm{g0}`, where :math:`\boldsymbol{u}_\mathrm{g0}` is the operating-point quantity. From :eq:`disturbance`--:eq:`frequency_tracking`, the linearized model for the estimation-error dynamics is obtained as
 
 .. math::
-    \frac{\mathrm{d}\Delta \tilde{\boldsymbol{u}}_\mathrm{g}}{\mathrm{d}t} &= -\alpha_\mathrm{g}\Delta \tilde{\boldsymbol{u}}_\mathrm{g} + \mathrm{j}\boldsymbol{u}_\mathrm{g0} (\Delta \omega_\mathrm{g} - \Delta \hat{\omega}_\mathrm{g}) \\
-    \frac{\mathrm{d}\Delta \hat{\omega}_\mathrm{g}}{\mathrm{d}t} &= k_\omega\mathrm{Im}\left\{ \frac{\Delta \tilde{\boldsymbol{u}}_\mathrm{g}}{\boldsymbol{u}_\mathrm{g0}} \right\} 
+    \frac{\mathrm{d}\Delta \tilde{\boldsymbol{u}}_\mathrm{g}}{\mathrm{d}t} &= -\alpha_\mathrm{o}\Delta \tilde{\boldsymbol{u}}_\mathrm{g} + \mathrm{j}\boldsymbol{u}_\mathrm{g0} (\Delta \omega_\mathrm{g} - \Delta \hat{\omega}_\mathrm{g}) \\
+    \frac{\mathrm{d}\Delta \hat{\omega}_\mathrm{g}}{\mathrm{d}t} &= k_\omega\mathrm{Im}\left\{ \frac{\Delta \tilde{\boldsymbol{u}}_\mathrm{g}}{\boldsymbol{u}_\mathrm{g0}} \right\}
     :label: linearized_model
 
-where :math:`\Delta \tilde{\boldsymbol{u}}_\mathrm{g} = \Delta\boldsymbol{u}_\mathrm{g} - \Delta \hat{\boldsymbol{u}}_\mathrm{g}` is the estimation error. 
+where :math:`\Delta \tilde{\boldsymbol{u}}_\mathrm{g} = \Delta\boldsymbol{u}_\mathrm{g} - \Delta \hat{\boldsymbol{u}}_\mathrm{g}` is the estimation error.
 
 First, assume that the grid frequency :math:`\omega_\mathrm{g}` is constant and the frequency tracking is disabled. From :eq:`linearized_model`, the closed-loop transfer function from the PCC voltage to its estimate becomes
 
 .. math::
-    \frac{\Delta\hat{\boldsymbol{u}}_\mathrm{g}(s)}{\Delta\boldsymbol{u}_\mathrm{g}(s)} = \frac{\alpha_\mathrm{g}}{s + \alpha_\mathrm{g}} 
+    \frac{\Delta\hat{\boldsymbol{u}}_\mathrm{g}(s)}{\Delta\boldsymbol{u}_\mathrm{g}(s)} = \frac{\alpha_\mathrm{o}}{s + \alpha_\mathrm{o}}
     :label: closed_loop_pll
 
-It can be realized that both the angle and magnitude of the PCC voltage estimate converge with the bandwidth :math:`\alpha_\mathrm{g}`. 
+It can be realized that both the angle and magnitude of the PCC voltage estimate converge with the bandwidth :math:`\alpha_\mathrm{o}`.
 
-Next, the frequency-tracking dynamics are also considered. From :eq:`linearized_model`, the closed-loop transfer function from the grid angular frequency to its estimate becomes 
+Next, the frequency-tracking dynamics are also considered. From :eq:`linearized_model`, the closed-loop transfer function from the grid angular frequency to its estimate becomes
 
 .. math::
-    \frac{\Delta\hat{\omega}_\mathrm{g}(s)}{\Delta\omega_\mathrm{g}(s)} 
-    = \frac{k_\omega}{s^2 + \alpha_\mathrm{g}s + k_\omega}
+    \frac{\Delta\hat{\omega}_\mathrm{g}(s)}{\Delta\omega_\mathrm{g}(s)}
+    = \frac{k_\omega}{s^2 + \alpha_\mathrm{o}s + k_\omega}
     :label: closed_loop_pll_frequency_tracking
 
-Choosing :math:`k_\omega = \alpha_\mathrm{pll}^2` and :math:`\alpha_\mathrm{g} = 2\alpha_\mathrm{pll}` yields the double pole at :math:`s = -\alpha_\mathrm{pll}`, where :math:`\alpha_\mathrm{pll}` is the frequency-tracking bandwidth.
+Choosing :math:`k_\omega = \alpha_\mathrm{pll}^2` and :math:`\alpha_\mathrm{o} = 2\alpha_\mathrm{pll}` yields the double pole at :math:`s = -\alpha_\mathrm{pll}`, where :math:`\alpha_\mathrm{pll}` is the frequency-tracking bandwidth.
 
-This PLL in estimated PCC coordinates is implemented in the class :class:`motulator.grid.control.PLL`. The :doc:`/grid_examples/grid_following/index` examples use the PLL to synchronize with the grid. 
+This PLL in estimated PCC coordinates is implemented in the class :class:`motulator.grid.control.PLL`. The :doc:`/grid_examples/grid_following/index` examples use the PLL to synchronize with the grid.
 
 .. rubric:: References
 
-.. [#Kau1997] Kaura and Blasko, "Operation of a phase locked loop system under distorted utility conditions," in IEEE Trans. Ind. Appl., 1997, https://doi.org/10.1109/28.567077
+.. [#Kau1997] Kaura, Blasko, "Operation of a phase locked loop system under distorted utility conditions," IEEE Trans. Ind. Appl., 1997, https://doi.org/10.1109/28.567077
 
 .. [#Fra1997] Franklin, Powell, Workman, "Digital Control of Dynamic Systems," 3rd ed., Menlo Park, CA: Addison-Wesley, 1997
 
 .. [#Nur2024] Nurminen, Mourouvin, Hinkkanen, Kukkola, "Multifunctional grid-forming converter control based on a disturbance observer, "IEEE Trans. Power Electron., 2024, https://doi.org/10.1109/TPEL.2024.3433503
 
-.. [#Har2021] Harnefors, Kukkola, Routimo, Hinkkanen, Wang, "A universal controller for grid-connected voltage-source converters," IEEE J. Emerg. Sel. Topics Power Electron., 2021, https://doi.org/10.1109/JESTPE.2020.3039407 
\ No newline at end of file
+.. [#Har2021] Harnefors, Kukkola, Routimo, Hinkkanen, Wang, "A universal controller for grid-connected voltage-source converters," IEEE J. Emerg. Sel. Topics Power Electron., 2021, https://doi.org/10.1109/JESTPE.2020.3039407
diff --git a/examples/drive/flux_vector/README.txt b/examples/drive/flux_vector/README.txt
index c1035fd7..b8780218 100644
--- a/examples/drive/flux_vector/README.txt
+++ b/examples/drive/flux_vector/README.txt
@@ -9,4 +9,4 @@ These examples demonstrate flux-vector control of electric machine drives [#Pel2
 
 .. [#Awa2019] Awan, Hinkkanen, Bojoi, Pellegrino, "Stator-flux-oriented control of synchronous motors: A systematic design procedure," IEEE Trans. Ind. Appl., 2019, https://doi.org/10.1109/TIA.2019.2927316
 
-.. [#Tii2024] Tiitinen, Hinkkanen, Harnefors, "Design framework for sensorless control of synchronous machine drives," IEEE Trans. Ind. Electron., 2024, https://doi.org/10.1109/TIE.2024.3429650 
+.. [#Tii2024] Tiitinen, Hinkkanen, Harnefors, "Design framework for sensorless control of synchronous machine drives," IEEE Trans. Ind. Electron., 2024, https://doi.org/10.1109/TIE.2024.3429650
diff --git a/examples/drive/vhz/plot_vhz_ctrl_6step_im_2kw.py b/examples/drive/vhz/plot_vhz_ctrl_6step_im_2kw.py
index d0cd096e..df6091ba 100644
--- a/examples/drive/vhz/plot_vhz_ctrl_6step_im_2kw.py
+++ b/examples/drive/vhz/plot_vhz_ctrl_6step_im_2kw.py
@@ -2,10 +2,10 @@
 2.2-kW induction motor, 6-step mode
 ===================================
 
-This example simulates V/Hz control of a 2.2-kW induction motor drive. The 
-six-step overmodulation is enabled, which increases the fundamental voltage as 
-well as the harmonics. Since the PWM is not synchronized with the stator 
-frequency, the harmonic content also depends on the ratio between the stator 
+This example simulates V/Hz control of a 2.2-kW induction motor drive. The
+six-step overmodulation is enabled, which increases the fundamental voltage as
+well as the harmonics. Since the PWM is not synchronized with the stator
+frequency, the harmonic content also depends on the ratio between the stator
 frequency and the sampling frequency.
 
 """
diff --git a/examples/drive/vhz/plot_vhz_ctrl_im_2kw.py b/examples/drive/vhz/plot_vhz_ctrl_im_2kw.py
index fbc8a15f..f248e4a2 100644
--- a/examples/drive/vhz/plot_vhz_ctrl_im_2kw.py
+++ b/examples/drive/vhz/plot_vhz_ctrl_im_2kw.py
@@ -3,7 +3,7 @@
 ====================================
 
 A diode bridge, stiff three-phase grid, and a DC link is modeled. The default
-parameters in this example yield open-loop V/Hz control. 
+parameters in this example yield open-loop V/Hz control.
 
 """
 # %%
diff --git a/examples/grid/grid_following/plot_gfl_10kva.py b/examples/grid/grid_following/plot_gfl_10kva.py
index 6800aa4e..e3816662 100644
--- a/examples/grid/grid_following/plot_gfl_10kva.py
+++ b/examples/grid/grid_following/plot_gfl_10kva.py
@@ -34,9 +34,7 @@
 
 # Create system model
 mdl = model.GridConverterSystem(converter, ac_filter, ac_source)
-
-# Uncomment line below to enable the PWM model
-# mdl.pwm = model.CarrierComparison()
+# mdl.pwm = model.CarrierComparison()  # Uncomment to enable the PWM model
 
 # %%
 # Configure the control system.
diff --git a/examples/grid/grid_forming/plot_gfm_obs_13kva.py b/examples/grid/grid_forming/plot_gfm_obs_13kva.py
index 9b680187..897e00cc 100644
--- a/examples/grid/grid_forming/plot_gfm_obs_13kva.py
+++ b/examples/grid/grid_forming/plot_gfm_obs_13kva.py
@@ -25,7 +25,7 @@
 
 # Filter and grid parameters
 par = ACFilterPars(L_fc=.15*base.L, R_fc=.05*base.Z, L_g=.74*base.L)
-# par.L_g = 0 # Uncomment this line to simulate a strong grid
+# par.L_g = 0  # Uncomment this line to simulate a strong grid
 ac_filter = model.ACFilter(par)
 ac_source = model.ThreePhaseVoltageSource(w_g=base.w, abs_e_g=base.u)
 # Inverter with constant DC voltage
@@ -44,8 +44,8 @@
     nom_u=base.u,
     nom_w=base.w,
     max_i=1.3*base.i,
-    T_s=100e-6,
-    R_a=.2*base.Z)
+    R_a=.2*base.Z,
+    T_s=100e-6)
 
 # Create the control system
 ctrl = control.ObserverBasedGridFormingControl(cfg)
diff --git a/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py b/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py
index 7f34e924..761d7ec0 100644
--- a/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py
+++ b/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py
@@ -22,7 +22,7 @@
 # Configure the system model.
 
 # Filter and grid
-par = ACFilterPars(L_fc=.15*base.L, L_g=.74*base.L)
+par = ACFilterPars(L_fc=.15*base.L, R_fc=.05*base.Z, L_g=.74*base.L)
 # par.L_g = 0  # Uncomment this line to simulate a strong grid
 ac_filter = model.ACFilter(par)
 # Grid voltage source with constant frequency and voltage magnitude
@@ -38,7 +38,12 @@
 
 # Control configuration parameters
 cfg = control.PowerSynchronizationControlCfg(
-    nom_u=base.u, nom_w=base.w, max_i=1.3*base.i, T_s=100e-6, R_a=.2*base.Z)
+    nom_u=base.u,
+    nom_w=base.w,
+    max_i=1.3*base.i,
+    R=.05*base.Z,
+    R_a=.2*base.Z,
+    T_s=100e-6)
 
 # Create the control system
 ctrl = control.PowerSynchronizationControl(cfg)
diff --git a/motulator/drive/control/im/_vhz.py b/motulator/drive/control/im/_vhz.py
index f3bc7c96..834a4ef3 100644
--- a/motulator/drive/control/im/_vhz.py
+++ b/motulator/drive/control/im/_vhz.py
@@ -42,7 +42,7 @@ class VHzControl(DriveControlSystem):
     """
     V/Hz control with the stator current feedback.
 
-    The method is similar to [#Hin2022]_. Open-loop V/Hz control can be 
+    The method is similar to [#Hin2022]_. Open-loop V/Hz control can be
     obtained as a special case by choosing::
 
         R_s, R_R = 0, 0
@@ -51,8 +51,8 @@ class VHzControl(DriveControlSystem):
     References
     ----------
     .. [#Hin2022] Hinkkanen, Tiitinen, Mölsä, Harnefors, "On the stability of
-       volts-per-hertz control for induction motors," IEEE J. Emerg. Sel. 
-       Topics Power Electron., 2022, 
+       volts-per-hertz control for induction motors," IEEE J. Emerg. Sel.
+       Topics Power Electron., 2022,
        https://doi.org/10.1109/JESTPE.2021.3060583
 
     """
diff --git a/motulator/drive/control/sm/_flux_vector.py b/motulator/drive/control/sm/_flux_vector.py
index c3a73fba..c7c5c6c4 100644
--- a/motulator/drive/control/sm/_flux_vector.py
+++ b/motulator/drive/control/sm/_flux_vector.py
@@ -15,11 +15,11 @@ class FluxVectorControl(DriveControlSystem):
     """
     Flux-vector control of synchronous machine drives.
 
-    This class implements a variant of flux-vector control [#Pel2009]_. Rotor 
-    coordinates as well as decoupling between the stator flux and torque 
-    channels are used according to [#Awa2019b]_. Here, the stator flux 
-    magnitude and the electromagnetic torque are selected as controllable 
-    variables, and proportional controllers are used for simplicity 
+    This class implements a variant of flux-vector control [#Pel2009]_. Rotor
+    coordinates as well as decoupling between the stator flux and torque
+    channels are used according to [#Awa2019b]_. Here, the stator flux
+    magnitude and the electromagnetic torque are selected as controllable
+    variables, and proportional controllers are used for simplicity
     [#Tii2024]_. The magnetic saturation is not considered in this
     implementation.
 
@@ -36,7 +36,7 @@ class FluxVectorControl(DriveControlSystem):
     alpha_o : float, optional
         Observer bandwidth (rad/s). The default is 2*pi*100.
     J : float, optional
-        Moment of inertia (kgm²). Needed only for the speed controller. 
+        Moment of inertia (kgm²). Needed only for the speed controller.
     T_s : float
         Sampling period (s). The default is 250e-6.
     sensorless : bool, optional
@@ -44,18 +44,18 @@ class FluxVectorControl(DriveControlSystem):
 
     References
     ----------
-    .. [#Pel2009] Pellegrino, Armando, Guglielmi, “Direct flux field-oriented 
-       control of IPM drives with variable DC link in the field-weakening 
-       region,” IEEE Trans. Ind. Appl., 2009, 
+    .. [#Pel2009] Pellegrino, Armando, Guglielmi, “Direct flux field-oriented
+       control of IPM drives with variable DC link in the field-weakening
+       region,” IEEE Trans. Ind. Appl., 2009,
        https://doi.org/10.1109/TIA.2009.2027167
 
-    .. [#Awa2019b] Awan, Hinkkanen, Bojoi, Pellegrino, "Stator-flux-oriented 
-       control of synchronous motors: A systematic design procedure," IEEE 
+    .. [#Awa2019b] Awan, Hinkkanen, Bojoi, Pellegrino, "Stator-flux-oriented
+       control of synchronous motors: A systematic design procedure," IEEE
        Trans. Ind. Appl., 2019, https://doi.org/10.1109/TIA.2019.2927316
 
-    .. [#Tii2024] Tiitinen, Hinkkanen, Harnefors, "Design framework for 
-       sensorless control of synchronous machine drives," IEEE Trans. Ind. 
-       Electron., 2024, https://doi.org/10.1109/TIE.2024.3429650 
+    .. [#Tii2024] Tiitinen, Hinkkanen, Harnefors, "Design framework for
+       sensorless control of synchronous machine drives," IEEE Trans. Ind.
+       Electron., 2024, https://doi.org/10.1109/TIE.2024.3429650
 
     """
 
@@ -127,7 +127,7 @@ class FluxTorqueReferenceCfg:
     par : SynchronousMachinePars
         Machine model parameters.
     max_i_s : float
-        Maximum stator current (A). 
+        Maximum stator current (A).
     min_psi_s : float, optional
         Minimum stator flux (Vs). The default is `par.psi_f`.
     max_psi_s : float, optional
diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index 50c4661a..6bf59d2d 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -14,18 +14,18 @@ class PLL:
     """
     Phase-locked loop including the voltage-magnitude filtering.
 
-    This class provides a simple frequency-tracking phase-locked loop. The 
-    magnitude of the measured PCC voltage is also filtered. 
+    This class provides a simple frequency-tracking phase-locked loop. The
+    magnitude of the measured PCC voltage is also filtered.
 
     Parameters
     ----------
     alpha_pll : float
         Frequency-tracking bandwidth (rad/s).
     abs_u_g0 : float
-        Initial value for the grid voltage estimate (V). 
+        Initial value for the grid voltage estimate (V).
     w_g0 : float
-        Initial value for the grid angular frequency estimate (rad/s). 
-        
+        Initial value for the grid angular frequency estimate (rad/s).
+
     """
 
     def __init__(self, alpha_pll, abs_u_g0, w_g0, theta_c0=0):
@@ -62,19 +62,19 @@ class DCBusVoltageController(PIController):
     """
     DC-bus voltage PI controller.
 
-    This controller regulates the energy stored in the DC-bus capacitor (scaled 
-    square of the DC-bus voltage) in order to have a linear closed-loop system 
-    [#Hur2001]_. 
+    This controller regulates the energy stored in the DC-bus capacitor (scaled
+    square of the DC-bus voltage) in order to have a linear closed-loop system
+    [#Hur2001]_.
 
     Parameters
     ----------
     C_dc : float
         DC-bus capacitance (F).
     alpha_dc : float
-        Approximate closed-loop bandwidth (rad/s). 
+        Approximate closed-loop bandwidth (rad/s).
     max_p : float, optional
         Limit for the maximum converter power (W). The default is `inf`.
-        
+
     References
     ----------
     .. [#Hur2001] Hur, Jung, Nam, "A fast dynamic DC-link power-balancing
@@ -100,10 +100,10 @@ def output(self, ref_u_dc, u_dc, u_ff=0):
 class GridConverterControlSystem(ControlSystem, ABC):
     """
     Base class for control of grid-connected converters.
-    
+
     This base class provides typical functionalities for control of
     grid-connected converters. This can be used both in power control and
-    DC-bus voltage control modes. 
+    DC-bus voltage control modes.
 
     Parameters
     ----------
@@ -117,7 +117,7 @@ class GridConverterControlSystem(ControlSystem, ABC):
 
             v : float | callable
                 Converter output voltage reference (V). Can be given either as
-                a constant or a function of time (s). 
+                a constant or a function of time (s).
             p_g : callable
                 Active power reference (W) as a function of time (s). This
                 signal is needed in power control mode.
@@ -142,7 +142,7 @@ def __init__(self, T_s):
     def get_electrical_measurements(self, fbk, mdl):
         """
         Measure the currents and voltages.
-        
+
         Parameters
         ----------
         fbk : SimpleNamespace
@@ -162,15 +162,11 @@ def get_electrical_measurements(self, fbk, mdl):
                 u_cs : complex
                     Realized converter output voltage (V) in stationary
                     coordinates. This signal is obtained from the PWM.
-                u_gs : complex
-                    PCC voltage (V) in stationary coordinates.
 
         """
         fbk.u_dc = mdl.converter.meas_dc_voltage()
         fbk.i_cs = abc2complex(mdl.ac_filter.meas_currents())
         fbk.u_cs = self.pwm.get_realized_voltage()
-        fbk.u_gs = abc2complex(mdl.ac_filter.meas_pcc_voltages())
-        # fbk.u_gs = abc2complex(mdl.ac_filter.meas_capacitor_voltages())
 
         return fbk
 
@@ -195,14 +191,14 @@ def get_power_reference(self, fbk, ref):
         Returns
         -------
         ref : SimpleNamespace
-            Reference signals, containing the following fields:
-                
+            Reference signals, containing some of the following fields:
+
                 u_dc : float
                     DC-bus voltage reference (V).
                 p_g : float
                     Active power reference (W).
                 q_g : float
-                    Reactive power reference (VAr).  
+                    Reactive power reference (VAr).
 
         """
         if self.dc_bus_voltage_ctrl:
@@ -211,12 +207,12 @@ def get_power_reference(self, fbk, ref):
             ref.p_g = self.dc_bus_voltage_ctrl.output(ref.u_dc, fbk.u_dc)
         else:
             # Power control mode
-            #ref.u_dc = None
             ref.p_g = self.ref.p_g(ref.t)
 
-        # Reactive power reference
-        ref.q_g = self.ref.q_g(ref.t) if callable(
-            self.ref.q_g) else self.ref.q_g
+        # Reactive power reference (if exists)
+        if hasattr(self.ref, 'q_g'):
+            ref.q_g = self.ref.q_g(ref.t) if callable(
+                self.ref.q_g) else self.ref.q_g
 
         return ref
 
diff --git a/motulator/grid/control/_grid_following.py b/motulator/grid/control/_grid_following.py
index 34f919c6..6ae1f8d6 100644
--- a/motulator/grid/control/_grid_following.py
+++ b/motulator/grid/control/_grid_following.py
@@ -7,6 +7,7 @@
 from motulator.common.control import ComplexPIController
 from motulator.grid.control._common import (
     CurrentLimiter, GridConverterControlSystem, PLL)
+from motulator.common.utils import abc2complex
 
 
 # %%
@@ -72,6 +73,8 @@ def __init__(self, cfg):
 
     def get_feedback_signals(self, mdl):
         fbk = super().get_feedback_signals(mdl)
+        fbk.u_gs = abc2complex(mdl.ac_filter.meas_pcc_voltages())
+        # fbk.u_gs = abc2complex(mdl.ac_filter.meas_capacitor_voltages())
         fbk = self.pll.output(fbk)
         s_g = 1.5*fbk.u_c*np.conj(fbk.i_c)
         fbk.p_g = s_g.real
diff --git a/motulator/grid/control/_observer_gfm.py b/motulator/grid/control/_observer_gfm.py
index 620093aa..ceff6707 100644
--- a/motulator/grid/control/_observer_gfm.py
+++ b/motulator/grid/control/_observer_gfm.py
@@ -5,7 +5,8 @@
 import numpy as np
 
 from motulator.common.utils import wrap
-from motulator.grid.control._common import GridConverterControlSystem
+from motulator.grid.control._common import (
+    GridConverterControlSystem, CurrentLimiter)
 
 
 # %%
@@ -28,12 +29,14 @@ class ObserverBasedGridFormingControlCfg:
         Total series resistance (Ω). The default is 0.
     R_a : float, optional
         Active resistance (Ω). The default is 0.25*nom_u/max_i.
-    T_s : float, optional
-        Sampling period (s). The default is 100e-6.
+    k_v : float, optional
+        Voltage control gain. The default is `alpha_o/nom_w`.
     alpha_c : float, optional
         Current control bandwidth (rad/s). The default is 2*pi*400.
     alpha_o : float, optional
         Observer gain (rad/s). The default is 2*pi*50.
+    T_s : float, optional
+        Sampling period (s). The default is 100e-6.
 
     """
     L: float
@@ -42,13 +45,16 @@ class ObserverBasedGridFormingControlCfg:
     max_i: float
     R: float = 0
     R_a: float = None
-    T_s: float = 100e-6
+    k_v: float = None
     alpha_c: float = 2*np.pi*400
     alpha_o: float = 2*np.pi*50
+    T_s: float = 100e-6
 
     def __post_init__(self):
         if self.R_a is None:
             self.R_a = .25*self.nom_u/self.max_i
+        if self.k_v is None:
+            self.k_v = self.alpha_o/self.nom_w
         self.k_c = self.alpha_c*self.L  # Current control gain
 
 
@@ -63,7 +69,7 @@ class ObserverBasedGridFormingControl(GridConverterControlSystem):
     Parameters
     ----------
     cfg : ObserverBasedGridFormingControlCfg
-        Controller configuration parameters.
+        Control system configuration parameters.
 
     Notes
     -----
@@ -83,9 +89,8 @@ class ObserverBasedGridFormingControl(GridConverterControlSystem):
     def __init__(self, cfg):
         super().__init__(cfg.T_s)
         self.cfg = cfg
-        self.observer = DisturbanceObserver(
-            cfg.nom_w, cfg.L, cfg.alpha_o, cfg.nom_u)
-        self.ref.q_g = 0
+        self.observer = DisturbanceObserver(cfg)
+        self.current_limiter = CurrentLimiter(cfg.max_i)
 
     def get_feedback_signals(self, mdl):
         """Get the feedback signals."""
@@ -108,21 +113,19 @@ def output(self, fbk):
             self.ref.v_c) else self.ref.v_c
 
         # Complex gains for grid-forming mode
-        abs_k_p = cfg.R_a/(1.5*ref.v_c)
-        abs_v_c = np.abs(fbk.v_c)
-        k_p = abs_k_p*fbk.v_c/abs_v_c if abs_v_c > 0 else 0
-        k_v = (1 - 1j)*fbk.v_c/abs_v_c if abs_v_c > 0 else 0
+        exp_j_theta = fbk.v_c/np.abs(fbk.v_c) if np.abs(fbk.v_c) > 0 else 1
+        k_p = cfg.R_a/(1.5*ref.v_c)*exp_j_theta
+        k_v = (1 - 1j*cfg.k_v)*exp_j_theta
 
         # Feedback correction for grid-forming mode
-        fbk.e_c = k_p*(ref.p_g - fbk.p_g) + k_v*(ref.v_c - abs_v_c)
+        fbk.e_c = k_p*(ref.p_g - fbk.p_g) + k_v*(ref.v_c - np.abs(fbk.v_c))
 
-        # Current limitation
+        # Transparent current limitation
         ref.i_c = fbk.i_c + fbk.e_c/cfg.k_c
-        if np.abs(ref.i_c) > cfg.max_i:
-            ref.i_c_lim = ref.i_c/np.abs(ref.i_c)*cfg.max_i
-            fbk.e_c = cfg.k_c*(ref.i_c_lim - fbk.i_c)
+        ref.i_c = self.current_limiter(ref.i_c)
+        fbk.e_c = cfg.k_c*(ref.i_c - fbk.i_c)
 
-        # Voltage reference
+        # Voltage reference (with optional resistive voltage drop compensation)
         ref.u_c = fbk.e_c + fbk.v_c + cfg.R*fbk.i_c
         ref.u_cs = np.exp(1j*fbk.theta_c)*ref.u_c
 
@@ -147,40 +150,30 @@ class DisturbanceObserver:
 
     Parameters
     ----------
-    w_g : float
-        Nominal grid angular frequency (rad/s).
-    L : float
-        Estimate of total inductance (H) between converter and grid.
-    alpha_o : float
-        Observer gain (rad/s).
-    v_c0 : float
-        Initial value of converter voltage state (V).
+    cfg : ObserverBasedGridFormingControlCfg
+        Control system configuration parameters.
 
     """
 
-    def __init__(self, w_g, L, alpha_o, v_c0):
-        self.w_g = w_g
-        self.w_c = w_g
-        self.L = L
-        self.k_o = alpha_o - 1j*w_g
-        # Initial states
-        self.v_cp = v_c0
+    def __init__(self, cfg):
+        self.alpha_o = cfg.alpha_o
+        self.L = cfg.L
+        self.w_g = cfg.nom_w
+        self.u_gp = cfg.nom_u
         self.theta_c = 0
 
     def output(self, fbk):
         """Compute the estimates."""
-        # Quasi-static converter voltage
-        fbk.v_c = self.v_cp - self.k_o*self.L*fbk.i_c
-        # Grid voltage
-        fbk.u_g = fbk.v_c - 1j*self.w_g*self.L*fbk.i_c
-        # Active and reactive power
-        fbk.p_g = 1.5*np.real(fbk.v_c*np.conj(fbk.i_c))
-        fbk.q_g = 1.5*np.imag(fbk.u_g*np.conj(fbk.i_c))
+        # Estimates for the quasi-static converter voltage and the grid voltage
+        fbk.v_c = self.u_gp - (self.alpha_o - 1j*self.w_g)*self.L*fbk.i_c
+        fbk.u_g = self.u_gp - self.alpha_o*self.L*fbk.i_c
+        # Active and reactive powers
+        s_g = 1.5*fbk.u_g*np.conj(fbk.i_c)
+        fbk.p_g, fbk.q_g = s_g.real, s_g.imag
         return fbk
 
     def update(self, fbk, ref):
-        """Update the observer states."""
-        self.v_cp += ref.T_s*(self.k_o + 1j*self.w_c)*fbk.e_c + 1j*(
-            self.w_g - self.w_c)*self.v_cp
-        self.theta_c += ref.T_s*self.w_c
+        """Update the states."""
+        self.u_gp += ref.T_s*self.alpha_o*fbk.e_c
+        self.theta_c += ref.T_s*self.w_g
         self.theta_c = wrap(self.theta_c)
diff --git a/motulator/grid/control/_power_synchronization.py b/motulator/grid/control/_power_synchronization.py
index 6a677b1e..fb943ee9 100644
--- a/motulator/grid/control/_power_synchronization.py
+++ b/motulator/grid/control/_power_synchronization.py
@@ -23,20 +23,23 @@ class PowerSynchronizationControlCfg:
         Nominal grid angular frequency (rad/s).
     max_i : float
         Maximum current (A), peak value.
+    R : float, optional
+        Total series resistance (Ω). The default is 0.
     R_a : float, optional
         Active resistance (Ω). The default is 0.25*nom_u/max_i.
-    T_s : float, optional
-        Sampling period (s). The default is 100e-6.
     w_b : float, optional
         Low-pass filter bandwidth (rad/s). The default is 2*pi*5.
+    T_s : float, optional
+        Sampling period (s). The default is 100e-6.
 
     """
     nom_u: float
     nom_w: float
     max_i: float
+    R: float = 0
     R_a: float = None
-    T_s: float = 100e-6
     w_b: float = 2*np.pi*5
+    T_s: float = 100e-6
 
     def __post_init__(self):
         if self.R_a is None:
@@ -60,8 +63,8 @@ class PowerSynchronizationControl(GridConverterControlSystem):
     References
     ----------
     .. [#Har2020] Harnefors, Rahman, Hinkkanen, Routimo, "Reference-feedforward
-        power-synchronization control," IEEE Trans. Power Electron., 2020,
-        https://doi.org/10.1109/TPEL.2020.2970991
+       power-synchronization control," IEEE Trans. Power Electron., 2020,
+       https://doi.org/10.1109/TPEL.2020.2970991
 
     """
 
@@ -69,7 +72,6 @@ def __init__(self, cfg):
         super().__init__(cfg.T_s)
         self.cfg = cfg
         self.current_limiter = CurrentLimiter(cfg.max_i)
-        self.ref.q_g = 0
         self.theta_c = 0
         self.i_c_flt = 0j
 
@@ -82,7 +84,8 @@ def get_feedback_signals(self, mdl):
         fbk.i_c = np.exp(-1j*fbk.theta_c)*fbk.i_cs
         fbk.u_c = np.exp(-1j*fbk.theta_c)*fbk.u_cs
         # Active power
-        fbk.p_c = 1.5*np.real(fbk.u_c*np.conj(fbk.i_c))
+        p_loss = 1.5*self.cfg.R*np.abs(fbk.i_c)**2
+        fbk.p_g = 1.5*np.real(fbk.u_c*np.conj(fbk.i_c)) - p_loss
         return fbk
 
     def output(self, fbk):
@@ -96,7 +99,7 @@ def output(self, fbk):
             self.ref.v_c) else self.ref.v_c
 
         # Calculation of power droop
-        fbk.w_c = cfg.nom_w + cfg.k_p_psc*(ref.p_g - fbk.p_c)
+        fbk.w_c = cfg.nom_w + cfg.k_p_psc*(ref.p_g - fbk.p_g)
 
         # Optionally, use of reference feedforward for d-axis current
         ref.i_c = ref.p_g/(1.5*ref.v_c) + 1j*fbk.i_c_flt.imag
@@ -106,7 +109,7 @@ def output(self, fbk):
         ref.i_c = self.current_limiter(ref.i_c)
 
         # Calculation of converter voltage output reference
-        ref.u_c = ref.v_c + cfg.R_a*(ref.i_c - fbk.i_c)
+        ref.u_c = ref.v_c + cfg.R_a*(ref.i_c - fbk.i_c) + cfg.R*fbk.i_c
         ref.u_cs = np.exp(1j*fbk.theta_c)*ref.u_c
 
         # Duty ratios for PWM
diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py
index 2b5d541b..c02d4d8d 100644
--- a/motulator/grid/utils/_plots.py
+++ b/motulator/grid/utils/_plots.py
@@ -182,12 +182,13 @@ def plot(sim, base=None, plot_pcc_voltage=True, plot_w=False, t_span=None):
         "--",
         label=r"$p_\mathrm{g,ref}$",
         ds="steps-post")
-    ax1.plot(
-        ctrl.t,
-        ctrl.ref.q_g/base.p,
-        "--",
-        label=r"$q_\mathrm{g,ref}$",
-        ds="steps-post")
+    if hasattr(ctrl.ref, "q_g"):
+        ax1.plot(
+            ctrl.t,
+            ctrl.ref.q_g/base.p,
+            "--",
+            label=r"$q_\mathrm{g,ref}$",
+            ds="steps-post")
     ax1.legend()
     ax1.set_xlim(t_span)
     ax1.set_xticklabels([])
diff --git a/pyproject.toml b/pyproject.toml
index 0c4b9719..3a94e826 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,6 +40,7 @@ dev = [
     "yapf",
     "pylint",
     "toml",
+    "pre-commit",
 ]
 doc = [
     "numpy",