From 4096ddcf88ff7bb6f1a93bccda2ef190f26094b0 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Fri, 7 Jun 2024 12:38:31 +0300
Subject: [PATCH 01/26] added initial grid directories

---
 ..._following_control_grid_converter_10kVA.py | 81 +++++++++++++++++++
 motulator/grid/__init__.py                    |  0
 motulator/grid/control/__init__.py            |  0
 motulator/grid/control/_grid_following.py     |  0
 motulator/grid/model/__init__.py              |  0
 motulator/grid/model/_const_freq_model.py     |  0
 6 files changed, 81 insertions(+)
 create mode 100644 examples/grid/plot_grid_following_control_grid_converter_10kVA.py
 create mode 100644 motulator/grid/__init__.py
 create mode 100644 motulator/grid/control/__init__.py
 create mode 100644 motulator/grid/control/_grid_following.py
 create mode 100644 motulator/grid/model/__init__.py
 create mode 100644 motulator/grid/model/_const_freq_model.py

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
new file mode 100644
index 000000000..b18e92c9d
--- /dev/null
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -0,0 +1,81 @@
+"""
+10-kVA grid following converter, power control
+==============================================
+    
+This example simulates a grid following controlled converter connected to a
+strong grid. The control system includes a phase-locked loop (PLL) to
+synchronize with the grid, a current reference generatior and a PI-based
+current controller.
+
+"""
+
+# %%
+# Imports.
+
+import numpy as np
+from motulator.grid import model, control
+from motulator import BaseValuesElectrical, plot_grid
+
+# To check the computation time of the program
+import time
+start_time = time.time()
+
+
+# %%
+# Compute base values based on the nominal values (just for figures).
+
+base_values = BaseValuesElectrical(
+    U_nom=400, I_nom=14.5, f_nom=50.0, P_nom=10e3)
+
+
+# %%
+# Create the system model.
+
+# grid impedance and filter model
+grid_filter = model.LFilter(L_f=10e-3, L_g=0, R_g=0)
+# AC grid model (either constant frequency or dynamic electromechanical model)
+grid_model = model.StiffSource(w_N=2*np.pi*50)
+converter = model.Inverter(u_dc=650)
+mdl = model.ac_grid.StiffSourceAndLFilterModel(
+    grid_filter, grid_model, converter)
+
+
+# %%
+# Configure the control system.
+
+# Control parameters
+pars = control.grid_following.GridFollowingCtrlPars(
+            L_f=10e-3,
+            f_sw = 5e3,
+            T_s = 1/(10e3),
+            i_max = 1.5*base_values.i,
+            )
+ctrl = control.grid_following.GridFollowingCtrl(pars)
+
+
+# %%
+# Set the time-dependent reference and disturbance signals.
+
+# Set the active and reactive power references
+ctrl.p_g_ref = lambda t: (t > .02)*(5e3)
+ctrl.q_g_ref = lambda t: (t > .04)*(4e3)
+
+# AC-voltage magnitude (to simulate voltage dips or short-circuits)
+e_g_abs_var =  lambda t: np.sqrt(2/3)*400
+mdl.grid_model.e_g_abs = e_g_abs_var # grid voltage magnitude
+
+
+# %%
+# Create the simulation object and simulate it.
+
+sim = model.Simulation(mdl, ctrl, pwm=False)
+sim.simulate(t_stop = .1)
+
+# Print the execution time
+print('\nExecution time: {:.2f} s'.format((time.time() - start_time)))
+
+
+# %%
+# Plot results in SI or per unit values.
+
+plot_grid(sim, base=base_values,plot_pcc_voltage=True)
diff --git a/motulator/grid/__init__.py b/motulator/grid/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/motulator/grid/control/__init__.py b/motulator/grid/control/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/motulator/grid/control/_grid_following.py b/motulator/grid/control/_grid_following.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/motulator/grid/model/__init__.py b/motulator/grid/model/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/motulator/grid/model/_const_freq_model.py b/motulator/grid/model/_const_freq_model.py
new file mode 100644
index 000000000..e69de29bb

From 3d0da54abdf709ea815d6f55678ba6e2ff6624fa Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Fri, 7 Jun 2024 14:01:41 +0300
Subject: [PATCH 02/26] grid/utils/ initializing, plots still have Bunch

---
 ..._following_control_grid_converter_10kVA.py |   7 +-
 motulator/grid/utils/__init__.py              |  12 +
 motulator/grid/utils/_helpers.py              |  98 +++++++
 motulator/grid/utils/_plots.py                | 248 ++++++++++++++++++
 4 files changed, 361 insertions(+), 4 deletions(-)
 create mode 100644 motulator/grid/utils/__init__.py
 create mode 100644 motulator/grid/utils/_helpers.py
 create mode 100644 motulator/grid/utils/_plots.py

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index b18e92c9d..8a7755d55 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -14,7 +14,7 @@
 
 import numpy as np
 from motulator.grid import model, control
-from motulator import BaseValuesElectrical, plot_grid
+from motulator.grid.utils import BaseValues, NominalValues, plot_grid
 
 # To check the computation time of the program
 import time
@@ -24,9 +24,8 @@
 # %%
 # Compute base values based on the nominal values (just for figures).
 
-base_values = BaseValuesElectrical(
-    U_nom=400, I_nom=14.5, f_nom=50.0, P_nom=10e3)
-
+nom = NominalValues(U=400, I=14.5, f=50, P=10e3)
+base = BaseValues.from_nominal(nom)
 
 # %%
 # Create the system model.
diff --git a/motulator/grid/utils/__init__.py b/motulator/grid/utils/__init__.py
new file mode 100644
index 000000000..5801eb5f7
--- /dev/null
+++ b/motulator/grid/utils/__init__.py
@@ -0,0 +1,12 @@
+"""This module contains utility functions for grid converters."""
+from motulator.grid.utils._helpers import BaseValues, NominalValues
+from motulator.grid.utils._plots import plot_grid
+from motulator.common.utils import Sequence, Step
+
+__all__ = [
+    "BaseValues",
+    "NominalValues",
+    "plot_grid",
+    "Sequence",
+    "Step",
+]
\ No newline at end of file
diff --git a/motulator/grid/utils/_helpers.py b/motulator/grid/utils/_helpers.py
new file mode 100644
index 000000000..bd44aa95e
--- /dev/null
+++ b/motulator/grid/utils/_helpers.py
@@ -0,0 +1,98 @@
+"""Helper functions and classes."""
+
+# %%
+from dataclasses import dataclass
+import numpy as np
+
+# %%
+@dataclass 
+class NominalValues:
+    """
+    Nominal values.
+
+    Parameters
+    ----------
+    U : float
+        Voltage (V, rms, line-line).
+    I : float
+        Current (A, rms).
+    f : float
+        Frequency (Hz).
+    P : float
+        Power (W).
+
+
+    """
+
+    U: float
+    I: float
+    f: float
+    P: float
+
+# %%
+@dataclass
+class BaseValues:
+    """
+    Base values.
+
+    Base values are computed from the nominal values and the number of pole
+    pairs. They can be used, e.g., for scaling the plotted waveforms.
+
+    Parameters
+    ----------
+    
+    u : float
+        Base voltage (V, peak, line-neutral).
+    i : float
+        Base current (A, peak).
+    w : float
+        Base angular frequency (rad/s).
+    psi : float
+        Base flux linkage (Vs).
+    p : float
+        Base power (W).
+    Z : float
+        Base impedance (Ω).
+    L : float
+        Base inductance (H).
+
+    """
+    # pylint: disable=too-many-instance-attributes
+    u: float
+    i: float
+    f: float
+    p: float
+
+@classmethod
+def from_nominal(cls, nom):
+     """
+        Compute base values from nominal values.
+
+        Parameters
+        ----------
+        nom : NominalValues
+            Nominal values containing the following fields:
+
+                U : float
+                    Voltage (V, rms, line-line).
+                I : float
+                    Current (A, rms).
+                f : float
+                    Frequency (Hz).
+
+
+        Returns
+        -------
+        BaseValues
+            Base values.
+
+        """
+    u = np.sqrt(2/3)*nom.U
+    i = np.sqrt(2)*nom.I
+    w = 2*np.pi*nom.f
+    psi = u/w
+    p = 1.5*u*i
+    Z = u/i
+    L = Z/w
+
+    return cls(u=u, i=i, w=w, psi=psi, p=p, Z=Z, L=L)
\ No newline at end of file
diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py
new file mode 100644
index 000000000..b4a41b82b
--- /dev/null
+++ b/motulator/grid/utils/_plots.py
@@ -0,0 +1,248 @@
+"""Example plotting scripts."""
+
+# %%
+import numpy as np
+import matplotlib.pyplot as plt
+from cycler import cycler
+
+from motulator.common.utils import complex2abc
+from gritulator._utils import Bunch
+
+# Plotting parameters
+plt.rcParams['axes.prop_cycle'] = cycler(color='brgcmyk')
+plt.rcParams['lines.linewidth'] = 1.
+plt.rcParams['axes.grid'] = True
+plt.rcParams.update({"text.usetex": False})
+
+
+# %%
+# pylint: disable=too-many-branches
+def plot_grid(
+        sim, t_range=None, base=None, plot_pcc_voltage=False, plot_w=False):
+    """
+    Plot example figures of grid converter simulations.
+
+    Plots figures in per-unit values, if the base values are given. Otherwise
+    SI units are used.
+
+    Parameters
+    ----------
+    sim : Simulation object
+        Should contain the simulated data.
+    t_range : 2-tuple, optional
+        Time range. The default is (0, sim.ctrl.t[-1]).
+    base : BaseValues, optional
+        Base values for scaling the waveforms.
+    plot_pcc_voltage : Boolean, optional
+        'True' if the user wants to plot the 3-phase waveform at the PCC. This
+        is an optional feature and the grid voltage is plotted by default.
+    plot_w : Boolean, optional
+        'True' if the user wants to plot the grid frequency instead of the
+        phase angles (by default).
+
+    """
+    FS = 16 # Font size of the plots axis
+    FL = 16 # Font size of the legends only
+    LW = 3 # Line width in plots
+    
+    
+    mdl = sim.mdl.data      # Continuous-time data
+    ctrl = sim.ctrl.data    # Discrete-time data
+    
+    # Check if the time span was given
+    if t_range is None:
+        t_range = (0, mdl.t[-1]) # Time span
+
+    # Check if the base values were given
+    if base is None:
+        pu_vals = False
+        # Scaling with unity base values except for power use kW
+        base = Bunch(w=1, u=1, i=1, p=1000)
+    else:
+        pu_vals = True
+        
+    # 3-phase quantities
+    i_g_abc = complex2abc(mdl.i_gs).T
+    u_g_abc = complex2abc(mdl.u_gs).T
+    e_g_abc = complex2abc(mdl.e_gs).T
+    
+    # grid voltage magnitude
+    abs_e_g = np.abs(mdl.e_gs)
+    
+    # Calculation of active and reactive powers
+    # p_g = 1.5*np.asarray(np.real(ctrl.u_g*np.conj(ctrl.i_c)))
+    # q_g = 1.5*np.asarray(np.imag(ctrl.u_g*np.conj(ctrl.i_c)))
+    p_g = 1.5*np.asarray(np.real(mdl.u_gs*np.conj(mdl.i_gs)))
+    q_g = 1.5*np.asarray(np.imag(mdl.u_gs*np.conj(mdl.i_gs)))
+    p_g_ref = np.asarray(ctrl.p_g_ref)
+    q_g_ref = np.asarray(ctrl.q_g_ref)
+    
+    # %%
+    fig, (ax1, ax2,ax3) = plt.subplots(3, 1, figsize=(10, 7))
+
+    if sim.ctrl.on_v_dc==False:
+        if plot_pcc_voltage == False:
+            # Subplot 1: Grid voltage
+            ax1.plot(mdl.t, e_g_abc/base.u, linewidth=LW)
+            ax1.legend([r'$e_g^a$',r'$e_g^b$',r'$e_g^c$'],
+                       prop={'size': FL}, loc= 'upper right')
+            ax1.set_xlim(t_range)
+            ax1.set_xticklabels([])
+        else:
+            # Subplot 1: PCC voltage
+            ax1.plot(mdl.t, u_g_abc/base.u, linewidth=LW)
+            ax1.legend([r'$u_g^a$',r'$u_g^b$',r'$u_g^c$'],
+                       prop={'size': FL}, loc= 'upper right')
+            ax1.set_xlim(t_range)
+            ax1.set_xticklabels([])
+            #ax1.set_ylabel('Grid voltage (V)')
+    else:
+        # Subplot 1: DC-bus voltage
+        ax1.plot(mdl.t, mdl.u_dc/base.u, linewidth=LW)
+        ax1.plot(ctrl.t, ctrl.u_dc_ref/base.u, '--', linewidth=LW)
+        ax1.legend([r'$u_{dc}$',r'$u_{dc,ref}$'],
+                   prop={'size': FL}, loc= 'upper right')
+        ax1.set_xlim(t_range)
+        ax1.set_xticklabels([])
+        #ax1.set_ylabel('DC-bus voltage (V)')
+    
+    # Subplot 2: Converter currents
+    ax2.plot(mdl.t, i_g_abc/base.i, linewidth=LW)
+    ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
+               ,prop={'size': FL}, loc= 'upper right')
+    ax2.set_xlim(t_range)
+    ax2.set_xticklabels([])
+    
+    if plot_w:
+        # Subplot 3: Grid and converter frequencies
+        ax3.plot(mdl.t, mdl.w_g/base.w, linewidth=LW)
+        ax3.plot(ctrl.t, ctrl.w_c/base.w, '--', linewidth=LW)
+        ax3.legend([r'$\omega_{g}$',r'$\omega_{c}$']
+                   ,prop={'size': FL}, loc= 'upper right')
+        ax3.set_xlim(t_range)
+    else:
+        # Subplot 3: Phase angles
+        ax3.plot(mdl.t, mdl.theta, linewidth=LW)
+        ax3.plot(ctrl.t, ctrl.theta_c, '--', linewidth=LW)
+        ax3.legend([r'$\theta_{g}$',r'$\theta_{c}$']
+                   ,prop={'size': FL}, loc= 'upper right')
+        ax3.set_xlim(t_range)
+
+    # Add axis labels
+    if pu_vals:
+        ax1.set_ylabel('Voltage (p.u.)')
+        ax2.set_ylabel('Current (p.u.)')
+    else:
+        ax1.set_ylabel('Voltage (V)')
+        ax2.set_ylabel('Current (A)')
+    if plot_w==False: 
+        ax3.set_ylabel('Angle (rad)')
+    elif plot_w and pu_vals:
+        ax3.set_ylabel('Frequency (p.u.)')
+    elif pu_vals == False:
+        ax3.set_ylabel('Frequency (rad/s)')
+    ax3.set_xlabel('Time (s)')
+
+    # Change font size
+    for item in ([ax1.title, ax1.xaxis.label, ax1.yaxis.label] +
+             ax1.get_xticklabels() + ax1.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    for item in ([ax2.title, ax2.xaxis.label, ax2.yaxis.label] +
+             ax2.get_xticklabels() + ax2.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    for item in ([ax3.title, ax3.xaxis.label, ax3.yaxis.label] +
+             ax3.get_xticklabels() + ax3.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    fig.align_ylabels()
+    plt.tight_layout()
+    plt.grid()
+    ax3.grid()
+    #plt.show()
+
+    
+    # %%
+    # Second figure
+    fig, (ax1, ax2,ax3) = plt.subplots(3, 1, figsize=(10, 7))
+
+    # Subplot 1: Active and reactive power
+    ax1.plot(mdl.t, p_g/base.p, linewidth=LW)
+    ax1.plot(mdl.t, q_g/base.p, linewidth=LW)
+    ax1.plot(ctrl.t, (p_g_ref/base.p), '--', linewidth=LW)
+    ax1.plot(ctrl.t, (q_g_ref/base.p), '--', linewidth=LW)
+    ax1.legend([r'$p_{g}$',r'$q_{g}$',r'$p_{g,ref}$',r'$q_{g,ref}$'],
+               prop={'size': FL}, loc= 'upper right')
+    ax1.set_xlim(t_range)
+    ax1.set_xticklabels([])
+    
+    # Subplot 2: Converter currents
+    ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
+    ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
+    ax2.plot(ctrl.t, np.real(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+    ax2.plot(ctrl.t, np.imag(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+    #ax2.plot(mdl.t, mdl.iL, linewidth=LW) converter-side dc current for debug
+    ax2.legend([r'$i_{c}^d$',r'$i_{c}^q$',r'$i_{c,ref}^d$',r'$i_{c,ref}^q$'],
+               prop={'size': FL}, loc= 'upper right')
+    ax2.set_xlim(t_range)
+    ax2.set_xticklabels([])
+
+    # Subplot 3: Converter voltage reference and grid voltage
+    ax3.plot(ctrl.t,np.real(ctrl.u_c/base.u), 
+            ctrl.t,np.imag(ctrl.u_c/base.u), linewidth=LW)
+    ax3.plot(mdl.t,np.real(abs_e_g/base.u),'k--', 
+             linewidth=LW)
+    ax3.legend([r'$u_{c}^d$', r'$u_{c}^q$', r'$|e_{g}|$'], 
+                prop={'size': FS}, loc= 'upper right')
+    ax3.set_xlim(t_range)
+    
+    # Change font size
+    for item in ([ax2.title, ax2.xaxis.label, ax2.yaxis.label] +
+             ax2.get_xticklabels() + ax2.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    for item in ([ax1.title, ax1.xaxis.label, ax1.yaxis.label] +
+             ax1.get_xticklabels() + ax1.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    for item in ([ax3.title, ax3.xaxis.label, ax3.yaxis.label] +
+             ax3.get_xticklabels() + ax3.get_yticklabels()):
+        item.set_fontsize(FS)
+
+    # Add axis labels
+    if pu_vals:
+        ax1.set_ylabel('Power (p.u.)')
+        ax2.set_ylabel('Current (p.u.)')
+        ax3.set_ylabel('Voltage (p.u.)')
+    else:
+        ax1.set_ylabel('Power (kW, kVar)')
+        ax2.set_ylabel('Current (A)')
+        ax3.set_ylabel('Voltage (V)')
+    ax3.set_xlabel('Time (s)')
+
+    fig.align_ylabels()
+    plt.tight_layout()
+    plt.grid()
+    ax3.grid()
+    plt.show()    
+
+
+
+# %%
+def save_plot(name):
+    """
+    Save figures.
+
+    This saves figures in a folder "figures" in the current directory. If the
+    folder doesn't exist, it is created.
+
+    Parameters
+    ----------
+    name : string
+        Name for the figure
+    plt : object
+        Handle for the figure to be saved
+
+    """
+    plt.savefig(name + '.pdf')

From 67f7e3cdb7a198e9df55e08558724b759421f5df Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Fri, 7 Jun 2024 14:35:05 +0300
Subject: [PATCH 03/26] grid/utils init & helpers done

---
 ..._following_control_grid_converter_10kVA.py |  2 +-
 motulator/grid/utils/__init__.py              |  2 +-
 motulator/grid/utils/_helpers.py              | 35 +++++++++++--------
 3 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 8a7755d55..a3b934188 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -77,4 +77,4 @@
 # %%
 # Plot results in SI or per unit values.
 
-plot_grid(sim, base=base_values,plot_pcc_voltage=True)
+plot_grid(sim, base, plot_pcc_voltage=True)
diff --git a/motulator/grid/utils/__init__.py b/motulator/grid/utils/__init__.py
index 5801eb5f7..5a76a4a5f 100644
--- a/motulator/grid/utils/__init__.py
+++ b/motulator/grid/utils/__init__.py
@@ -9,4 +9,4 @@
     "plot_grid",
     "Sequence",
     "Step",
-]
\ No newline at end of file
+]
diff --git a/motulator/grid/utils/_helpers.py b/motulator/grid/utils/_helpers.py
index bd44aa95e..f7db47a21 100644
--- a/motulator/grid/utils/_helpers.py
+++ b/motulator/grid/utils/_helpers.py
@@ -5,7 +5,7 @@
 import numpy as np
 
 # %%
-@dataclass 
+@dataclass
 class NominalValues:
     """
     Nominal values.
@@ -49,6 +49,8 @@ class BaseValues:
         Base angular frequency (rad/s).
     psi : float
         Base flux linkage (Vs).
+    f : float
+        Base frequency (Hz).
     p : float
         Base power (W).
     Z : float
@@ -60,12 +62,17 @@ class BaseValues:
     # pylint: disable=too-many-instance-attributes
     u: float
     i: float
+    w: float
+    psi: float
     f: float
     p: float
+    f: float
+    Z: float
+    L: float
 
-@classmethod
-def from_nominal(cls, nom):
-     """
+    @classmethod
+    def from_nominal(cls, nom):
+        """
         Compute base values from nominal values.
 
         Parameters
@@ -84,15 +91,13 @@ def from_nominal(cls, nom):
         Returns
         -------
         BaseValues
-            Base values.
-
         """
-    u = np.sqrt(2/3)*nom.U
-    i = np.sqrt(2)*nom.I
-    w = 2*np.pi*nom.f
-    psi = u/w
-    p = 1.5*u*i
-    Z = u/i
-    L = Z/w
-
-    return cls(u=u, i=i, w=w, psi=psi, p=p, Z=Z, L=L)
\ No newline at end of file
+        u = np.sqrt(2/3)*nom.U
+        i = np.sqrt(2)*nom.I
+        w = 2*np.pi*nom.f
+        psi = u/w
+        p = 1.5*u*i
+        Z = u/i
+        L = Z/w
+
+        return cls(u=u, i=i, w=w, psi=psi, p=p, Z=Z, L=L)

From b6e33fd20c83729fe3799f15ce9fbbe7c69a387f Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Fri, 7 Jun 2024 14:48:04 +0300
Subject: [PATCH 04/26] add gfl models from gritulator

---
 ..._following_control_grid_converter_10kVA.py |   2 +-
 motulator/common/control/__init__.py          |   4 +-
 motulator/common/control/_control.py          | 106 +++++
 motulator/grid/__init__.py                    |   1 +
 motulator/grid/control/__init__.py            |   2 +
 motulator/grid/control/_grid_following.py     | 404 ++++++++++++++++++
 motulator/grid/model/__init__.py              |  11 +
 motulator/grid/model/_const_freq_model.py     | 271 ++++++++++++
 motulator/grid/model/_grid_filter.py          | 298 +++++++++++++
 motulator/grid/model/_grid_volt_source.py     | 251 +++++++++++
 10 files changed, 1348 insertions(+), 2 deletions(-)
 create mode 100644 motulator/grid/model/_grid_filter.py
 create mode 100644 motulator/grid/model/_grid_volt_source.py

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index b18e92c9d..9526ceb5e 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -14,7 +14,7 @@
 
 import numpy as np
 from motulator.grid import model, control
-from motulator import BaseValuesElectrical, plot_grid
+from motulator.grid.utils import BaseValuesElectrical, plot_grid
 
 # To check the computation time of the program
 import time
diff --git a/motulator/common/control/__init__.py b/motulator/common/control/__init__.py
index 5cd391e7e..14ab48f7d 100644
--- a/motulator/common/control/__init__.py
+++ b/motulator/common/control/__init__.py
@@ -1,11 +1,13 @@
 """Common control functions and classes."""
 from motulator.common.control._control import (
-    Ctrl, ComplexPICtrl, PICtrl, PWM, RateLimiter)
+    Ctrl, ComplexPICtrl, ComplexFFPICtrl, PICtrl, PWM, RateLimiter, Clock, )
 
 __all__ = [
     "Ctrl",
     "ComplexPICtrl",
+    "ComplexFFPICtrl",
     "PICtrl",
     "PWM",
     "RateLimiter",
+    "Clock",
 ]
diff --git a/motulator/common/control/_control.py b/motulator/common/control/_control.py
index 35219507a..642809f7f 100644
--- a/motulator/common/control/_control.py
+++ b/motulator/common/control/_control.py
@@ -380,6 +380,112 @@ def update(self, T_s, u, w):
         self.u_i += T_s*(self.alpha_i + 1j*w)*(u - self.v)
 
 
+# %%
+class ComplexFFPICtrl:
+    """
+    2DOF Synchronous-frame complex-vector PI controller with feedforward.
+
+    This implements a discrete-time 2DOF synchronous-frame complex-vector PI
+    controller similar to [#Bri2000]_, with an additional feedforward signal.                        
+    The gain selection corresponding to internal-model-control (IMC) is
+    equivalent to the continuous-time version given in [#Har2009]_::
+
+        u = k_p*(i_ref - i) + k_i/s*(i_ref - i) + 1j*w*L_f*i + u_g_ff 
+                
+    where `u` is the controller output, `i_ref` is the reference signal, `i` is
+    the feedback signal, u_g_ff is the (filtered) feedforward signal, `w` is
+    the angular speed of synchronous coordinates, '1j*w*L_f' is the decoupling
+    term estimate, and `1/s` refers to integration. This algorithm is
+    compatible with both real and complex signals. The integrator anti-windup
+    is implemented based on the realized controller output.
+
+    Parameters
+    ----------
+    k_p : float
+        Proportional gain.
+    k_i : float
+        Integral gain.
+    k_t : float, optional
+        Reference-feedforward gain. The default is `k_p`.
+    L_f : float, optional
+        Synchronous frame decoupling gain. The default is 0.
+
+    Attributes
+    ----------
+    v : complex
+        Input disturbance estimate.
+    u_i : complex
+        Integral state.
+
+    Notes
+    -----
+    This contoller can be used, e.g., as a current controller. In this case,
+    `i` corresponds to the converter current and `u` to the converter voltage.
+    
+    References
+    ----------
+       
+    .. [#Har2009] Harnefors, Bongiorno, "Current controller design
+       for passivity of the input admittance," 2009 13th European Conference
+       on Power Electronics and Applications, Barcelona, Spain, 2009.
+
+    """
+
+    def __init__(self, k_p, k_i, k_t=None, L_f=None):
+        # Gains
+        self.k_p = k_p
+        self.k_t = k_t if k_t is not None else k_p
+        self.alpha_i = k_i/self.k_t  # Inverse of the integration time T_i
+        self.L_f = L_f if L_f is not None else 0
+        # States
+        self.v, self.u_i = 0, 0
+
+    def output(self, i_ref, i, u_ff, w):
+        """
+        Compute the controller output.
+
+        Parameters
+        ----------
+        i_ref : complex
+            Reference signal.
+        i : complex
+            Feedback signal.
+        u_ff : complex
+            Feedforward signal.
+        w : float
+            Angular speed of the reference frame (rad/s). 
+
+        Returns
+        -------
+        u : complex
+            Controller output.
+
+        """
+        # Disturbance input estimate
+        self.v = self.u_i - (self.k_p - self.k_t)*i + (u_ff
+        + 1j*w*self.L_f*i)
+
+        # Controller output
+        u = self.k_t*(i_ref - i) + self.v
+
+        return u
+
+    def update(self, T_s, u_lim):
+        """
+        Update the integral state.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+        u_lim : complex
+            Realized (limited) controller output. If the actuator does not
+            saturate, ``u_lim = u``.
+
+        """
+        self.u_i += T_s*self.alpha_i*(u_lim - self.v)
+
+
 # %%
 class RateLimiter:
     """
diff --git a/motulator/grid/__init__.py b/motulator/grid/__init__.py
index e69de29bb..3f890cb81 100644
--- a/motulator/grid/__init__.py
+++ b/motulator/grid/__init__.py
@@ -0,0 +1 @@
+"""Grid converter models and controls."""
diff --git a/motulator/grid/control/__init__.py b/motulator/grid/control/__init__.py
index e69de29bb..18bebfd8d 100644
--- a/motulator/grid/control/__init__.py
+++ b/motulator/grid/control/__init__.py
@@ -0,0 +1,2 @@
+"""Controls for grid-connected converters."""
+
diff --git a/motulator/grid/control/_grid_following.py b/motulator/grid/control/_grid_following.py
index e69de29bb..485e2718e 100644
--- a/motulator/grid/control/_grid_following.py
+++ b/motulator/grid/control/_grid_following.py
@@ -0,0 +1,404 @@
+"""grid following control methods for grid onverters."""
+
+# %%
+from __future__ import annotations
+from collections.abc import Callable
+from dataclasses import dataclass, field
+
+from types import SimpleNamespace
+#from motulator._utils import Bunch
+import numpy as np
+
+from motulator.common.utils._utils import abc2complex
+from motulator.common.control import (Ctrl, PWM, ComplexFFPICtrl, Clock)
+
+
+# %%
+@dataclass
+class GridFollowingCtrlPars:
+    """
+    grid-following control parameters.
+
+    """
+    # pylint: disable=too-many-instance-attributes
+    # General control parameters
+    p_g_ref: Callable[[float], float] = field(
+        repr=False, default=lambda t: 0) # active power reference
+    q_g_ref: Callable[[float], float] = field(
+        repr=False, default=lambda t: 0) # reactive power reference
+    u_dc_ref: Callable[[float], float] = field(
+        repr=False, default=lambda t: 650) # DC voltage reference, only used if
+                                    # the dc voltage control mode is activated.
+    T_s: float = 1/(16e3) # sampling time of the controller.
+    delay: int = 1
+    u_gN: float = np.sqrt(2/3)*400  # PCC voltage, in volts.
+    w_g: float = 2*np.pi*50 # grid frequency, in Hz
+    f_sw: float = 8e3 # switching frequency, in Hz.
+
+    # Scaling ratio of the abc/dq transformation
+    k_scal: float = 3/2
+
+    # Current controller parameters
+    alpha_c: float = 2*np.pi*400 # current controller bandwidth.
+
+    # Phase Locked Loop (PLL) control parameters
+    w_0_pll: float = 2*np.pi*20 # undamped natural frequency
+    zeta: float = 1 # damping ratio
+
+    # Low pass filter for voltage feedforward term
+    alpha_ff: float = 2*np.pi*(4*50) # low pass filter bandwidth
+
+    # Use the filter capacitance voltage measurement or PCC voltage
+    on_u_cap: bool = 0 # 1 if capacitor voltage is used, 0 if PCC is used
+
+    # DC-voltage controller
+    on_v_dc: bool = 0 # put 1 to activate dc voltage controller. 0 is p-mode
+    zeta_dc: float = 1 # damping ratio
+    p_max: float = 10e3 # maximum power reference, in W.
+    w_0_dc: float = 2*np.pi*30 # controller undamped natural frequency, in rad/s.
+
+    # Current limitation
+    i_max: float = 20 # maximum current modulus in A
+
+    # Passive component parameter estimates
+    L_f: float = 10e-3 # filter inductance, in H.
+    C_dc: float = 1e-3 # DC bus capacitance, in F.
+
+
+# %%
+class GridFollowingCtrl(Ctrl):
+    """
+    Grid following control for power converters.
+
+    Parameters
+    ----------
+    pars : GridFollowingCtrlPars
+        Control parameters.
+
+    """
+
+    # pylint: disable=too-many-instance-attributes
+    def __init__(self, pars):
+        super().__init__()
+        self.t = 0
+        self.T_s = pars.T_s
+        # Instantiate classes
+        self.pwm = PWM(six_step=False)
+        self.clock = Clock()
+        self.pll = PLL(pars)
+        self.current_ctrl = CurrentCtrl(pars, pars.alpha_c)
+        self.current_ref_calc = CurrentRefCalc(pars)
+        self.dc_bus_volt_ctrl = DCBusVoltCtrl(
+            pars.zeta_dc, pars.w_0_dc, pars.p_max)
+        # Parameters
+        self.u_gN = pars.u_gN
+        self.w_g = pars.w_g
+        self.f_sw = pars.f_sw
+        self.L_f = pars.L_f
+        # Power references
+        self.p_g_ref = pars.p_g_ref
+        self.q_g_ref = pars.q_g_ref
+        # Activation/deactivation of the DC voltage controller
+        self.on_v_dc = pars.on_v_dc
+        # DC voltage reference
+        self.u_dc_ref = pars.u_dc_ref
+        # DC-bus capacitance
+        self.C_dc = pars.C_dc
+        # Calculated current controller gains:
+        self.k_p_i = pars.alpha_c*pars.L_f
+        self.k_i_i = np.power(pars.alpha_c,2)*pars.L_f
+        self.r_i = pars.alpha_c*pars.L_f
+        # Calculated maximum current in A
+        self.i_max = pars.i_max
+        # Calculated PLL estimator gains
+        self.k_p_pll = 2*pars.zeta*pars.w_0_pll/pars.u_gN
+        self.k_i_pll = pars.w_0_pll*pars.w_0_pll/pars.u_gN
+        # Low pass filter for voltage feedforward term
+        self.alpha_ff = pars.alpha_ff
+        # Measure the voltage at the PCC or the capacitor of LCL (if used)?
+        self.on_u_cap = pars.on_u_cap
+        # States
+        self.u_c_i = 0j
+        self.theta_p = 0
+        self.u_g_filt = pars.u_gN + 1j*0
+
+    def __call__(self, mdl):
+        """
+        Run the main control loop.
+
+        Parameters
+        ----------
+        mdl : StiffSourceAndLFilterModel / StiffSourceAndLCLFilterModel /
+              ACFlexSourceAndLFilterModel / ACFlexSourceAndLCLFilterModel /
+              DCBusAndLFilterModel / DCBusAndLCLFilterModel
+            Continuous-time model of a voltage-source grid model with a filter
+            and a resistive-inductive line for getting the feedback signals.
+
+        Returns
+        -------
+        T_s : float
+            Sampling period.
+        d_abc_ref : ndarray, shape (3,)
+            Duty ratio references.
+
+        """
+        # Measure the feedback signals
+        i_c_abc = mdl.grid_filter.meas_currents()
+        u_dc = mdl.converter.meas_dc_voltage()
+        if self.on_u_cap == True:
+            u_g_abc = mdl.grid_filter.meas_cap_voltage()
+        else:
+            u_g_abc = mdl.grid_filter.meas_pcc_voltage()
+        # Obtain the converter voltage calculated with the PWM
+        u_c = self.pwm.realized_voltage
+
+        # Define the active and reactive power references at the given time 
+        u_dc_ref = self.u_dc_ref(self.clock.t)
+        # Definition of capacitance energy variables for the DC-bus controller
+        W_dc_ref = 0.5*self.C_dc*u_dc_ref**2
+        W_dc = 0.5*self.C_dc*u_dc**2
+        if self.on_v_dc:
+            p_g_ref = self.dc_bus_volt_ctrl.output(W_dc_ref, W_dc)
+            q_g_ref = self.q_g_ref(self.clock.t)
+        else:
+            p_g_ref = self.p_g_ref(self.clock.t)
+            q_g_ref = self.q_g_ref(self.clock.t)
+
+        # Generate the current references
+        i_c_ref = self.current_ref_calc.output(p_g_ref, q_g_ref)
+
+        # Transform the measured current in dq frame
+        i_c = np.exp(-1j*self.theta_p)*abc2complex(i_c_abc)
+
+        # Calculation of PCC voltage in synchronous frame
+        u_g = np.exp(-1j*self.theta_p)*abc2complex(u_g_abc)
+
+        # Use of PLL to bring ugq to zero
+        u_g_q, abs_u_g, w_pll, theta_pll = self.pll.output(u_g_abc)
+
+        # Calculation of the modulus of current reference
+        i_abs = np.abs(i_c_ref)
+        i_cd_ref = np.real(i_c_ref)
+        i_cq_ref = np.imag(i_c_ref)
+
+        # And current limitation algorithm
+        if i_abs > 0:
+            i_ratio = self.i_max/i_abs
+            i_cd_ref = np.sign(i_cd_ref)*np.min(
+                [i_ratio*np.abs(i_cd_ref),np.abs(i_cd_ref)])
+            i_cq_ref = np.sign(i_cq_ref)*np.min(
+                [i_ratio*np.abs(i_cq_ref),np.abs(i_cq_ref)])
+            i_c_ref = i_cd_ref + 1j*i_cq_ref
+
+        # Low pass filter for the feedforward PCC voltage:
+        u_g_filt = self.u_g_filt
+
+        # Voltage reference generation in synchronous coordinates
+        u_c_ref = self.current_ctrl.output(i_c_ref, i_c, u_g_filt, self.w_g)
+
+        # Use the function from control commons:
+        d_abc_ref = self.pwm(self.T_s, u_c_ref, u_dc,
+                                           self.theta_p, self.w_g)
+
+        # Data logging
+        data = SimpleNamespace(
+            w_c = w_pll, theta_c = self.theta_p, u_c_ref = u_c_ref,
+            u_c = u_c, i_c = i_c, abs_u_g = abs_u_g,
+            d_abc_ref = d_abc_ref, i_c_ref = i_c_ref, u_dc = u_dc,
+            t = self.clock.t, p_g_ref = p_g_ref, u_dc_ref = u_dc_ref,
+            q_g_ref=q_g_ref, u_g = u_g,
+                     )
+        self.save(data)
+
+        # Update the states
+        self.theta_p = theta_pll
+        self.w_p = w_pll
+        self.current_ctrl.update(self.T_s, u_c)
+        self.clock.update(self.T_s)
+        self.pll.update(u_g_q)
+        if self.on_v_dc:
+            self.dc_bus_volt_ctrl.update(self.T_s, p_g_ref)
+        # Update the low pass filer integrator for feedforward action
+        self.u_g_filt = (1 - self.T_s*self.alpha_ff)*u_g_filt + (
+            self.T_s*self.alpha_ff*u_g)
+
+        return self.T_s, d_abc_ref
+
+
+# %%
+class PLL:
+
+    """
+    PLL synchronizing loop.
+
+    Parameters
+    ----------
+    u_g_abc : ndarray, shape (3,)
+        Phase voltages at the PCC.
+
+    Returns
+    -------
+    u_g_q : float
+        q-axis of the PCC voltage (V)
+    abs_u_g : float
+        amplitude of the voltage waveform, in V
+    theta_pll : float
+        estimated phase angle (in rad).
+        
+    """
+
+    def __init__(self, pars):
+
+        """
+        Parameters
+        ----------
+        pars : GridFollowingCtrlPars
+           Control parameters.
+    
+        """
+        self.T_s = pars.T_s
+        self.w_0_pll = pars.w_0_pll
+        self.k_p_pll = 2*pars.zeta*pars.w_0_pll/pars.u_gN
+        self.k_i_pll = pars.w_0_pll*pars.w_0_pll/pars.u_gN
+        # Initial states
+        self.w_pll = pars.w_g
+        self.theta_p = 0
+
+
+    def output(self, u_g_abc):
+
+        """
+        Compute the estimated frequency and phase angle using the PLL.
+    
+        Parameters
+        ----------
+        u_g_abc : ndarray, shape (3,)
+            Grid 3-phase voltage.
+    
+        Returns
+        -------
+        u_g_q : float
+            Error signal (in V, corresponds to the q-axis grid voltage).
+        abs_u_g : float
+            magnitude of the grid voltage vector (in V).
+        w_g_pll : float
+            estimated grid frequency (in rad/s).
+        theta_pll : float
+            estimated phase angle (in rad).
+        """
+
+        u_g_ab = u_g_abc[0] - u_g_abc[1] # calculation of phase-to-phase voltages
+        u_g_bc = u_g_abc[1] - u_g_abc[2] # calculation of phase-to-phase voltages
+
+        # Calculation of u_g in complex form (stationary coordinates)
+        u_g_s = (2/3)*u_g_ab +(1/3)*u_g_bc + 1j*(np.sqrt(3)/(3))*u_g_bc
+        # And then in general coordinates
+        u_g = u_g_s*np.exp(-1j*self.theta_p)
+        # Definition of the error using the q-axis voltage
+        u_g_q = np.imag(u_g)
+
+        # Absolute value of the grid-voltage vector
+        abs_u_g = np.abs(u_g)
+
+        # Calculation of the estimated PLL frequency
+        w_g_pll = self.k_p_pll*u_g_q + self.w_pll
+
+        # Estimated phase angle
+        theta_pll = self.theta_p + self.T_s*w_g_pll
+
+        return u_g_q, abs_u_g, w_g_pll, theta_pll
+
+
+    def update(self, u_g_q):
+        """
+        Update the integral state.
+    
+        Parameters
+        ----------
+        u_g_q : real
+            Error signal (in V, corresponds to the q-axis grid voltage).
+    
+        """
+
+        # Calculation of the estimated PLL frequency
+        w_g_pll = self.k_p_pll*u_g_q + self.w_pll
+
+        # Update the integrator state
+        self.w_pll = self.w_pll + self.T_s*self.k_i_pll*u_g_q
+        # Update the grid-voltage angle state
+        self.theta_p = self.theta_p + self.T_s*w_g_pll
+        self.theta_p = np.mod(self.theta_p, 2*np.pi)    # Limit to [0, 2*pi]
+
+
+# %%
+class CurrentCtrl(ComplexFFPICtrl):
+    """
+    2DOF PI current controller for grid converters.
+
+    This class provides an interface for a current controller for grid 
+    converters. The gains are initialized based on the desired closed-loop 
+    bandwidth and the filter inductance. 
+
+    Parameters
+    ----------
+    par : ModelPars
+        Grid converter parameters, contains the filter inductance `L_f` (H).  
+    alpha_c : float
+        Closed-loop bandwidth (rad/s).
+
+    """
+
+    def __init__(self, par, alpha_c):
+        k_t = alpha_c*par.L_f
+        k_i = alpha_c*k_t
+        k_p = 2*k_t
+        L_f = par.L_f
+        super().__init__(k_p, k_i, k_t, L_f)
+
+
+# %%
+class CurrentRefCalc:
+
+    """
+    Current controller reference generator
+    
+    This class is used to generate the current references for the current
+    controllers based on the active and reactive power references.
+    
+    """
+
+    def __init__(self, pars):
+
+        """
+        Parameters
+        ----------
+        pars : GridFollowingCtrlPars
+            Control parameters.
+    
+        """
+        self.u_gN = pars.u_gN
+
+
+    def output(self, p_g_ref, q_g_ref):
+
+        """
+        Current reference genetator.
+    
+        Parameters
+        ----------
+        p_g_ref : float
+            active power reference
+        q_g_ref : float
+            reactive power reference
+    
+        Returns
+        -------
+        i_c_ref : float
+            current reference in the rotationary frame
+            
+        """
+
+        # Calculation of the current references in the stationary frame:
+        i_c_ref = 2*p_g_ref/(3*self.u_gN) -2*1j*q_g_ref/(3*self.u_gN)
+
+        return i_c_ref
diff --git a/motulator/grid/model/__init__.py b/motulator/grid/model/__init__.py
index e69de29bb..fd1e594ec 100644
--- a/motulator/grid/model/__init__.py
+++ b/motulator/grid/model/__init__.py
@@ -0,0 +1,11 @@
+"""Continuous-time grid converter interconnector models."""
+
+from motulator.grid.model._const_freq_model import (
+    StiffSourceAndLFilterModel,
+    StiffSourceAndLCLFilterModel,
+)
+
+__all__ = [
+    "StiffSourceAndLFilterModel",
+    "StiffSourceAndLCLFilterModel",
+]
diff --git a/motulator/grid/model/_const_freq_model.py b/motulator/grid/model/_const_freq_model.py
index e69de29bb..eca3d18c6 100644
--- a/motulator/grid/model/_const_freq_model.py
+++ b/motulator/grid/model/_const_freq_model.py
@@ -0,0 +1,271 @@
+# pylint: disable=C0103
+"""
+Electromechanical stiff AC grid and converter model interconnectors.
+
+    This interconnects the subsystems of a converter with a grid and provides
+    an interface to the solver. More complicated systems could be modeled using
+    a similar template. Peak-valued complex space vectors are used. The space
+    vector models are implemented in stationary coordinates.
+
+"""
+
+import numpy as np
+
+#from gritulator._utils import Bunch
+
+# %%
+class StiffSourceAndLFilterModel:
+    """
+    Continuous-time model for a stiff grid model with an RL impedance model.
+
+    Parameters
+    ----------
+    grid_filter : LFilter
+        RL line dynamic model.
+    grid_model : StiffSource
+        Constant voltage source model.
+    converter : Inverter | PWMInverter
+        Inverter model.
+
+    """
+
+    def __init__(self, grid_filter=None, grid_model=None, converter=None):
+        self.grid_filter = grid_filter
+        self.grid_model = grid_model
+        self.converter = converter
+
+        # Initial time
+        self.t0 = 0
+
+        # Store the solution in these lists
+        self.data = Bunch()
+        self.data.t, self.data.q = [], []
+        self.data.i_gs = []
+
+
+    def get_initial_values(self):
+        """
+        Get the initial values.
+
+        Returns
+        -------
+        x0 : complex list, length 1
+            Initial values of the state variables.
+
+        """
+        x0 = [self.grid_filter.i_gs0]
+
+        return x0
+
+    def set_initial_values(self, t0, x0):
+        """
+        Set the initial values.
+
+        Parameters
+        ----------
+        x0 : complex ndarray
+            Initial values of the state variables.
+
+        """
+        self.t0 = t0
+        self.grid_filter.i_gs0 = x0[0]
+        # calculation of converter-side voltage
+        u_cs0 = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
+        # calculation of grid-side voltage
+        e_gs0 = self.grid_model.voltages(t0)
+        # update pcc voltage
+        self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(x0[0], u_cs0, e_gs0)
+
+    def f(self, t, x):
+        """
+        Compute the complete state derivative list for the solver.
+
+        Parameters
+        ----------
+        t : float
+            Time.
+        x : complex ndarray
+            State vector.
+
+        Returns
+        -------
+        complex list
+            State derivatives.
+
+        """
+        # Unpack the states
+        i_gs = x
+        # Interconnections: outputs for computing the state derivatives
+        u_cs = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
+        e_gs = self.grid_model.voltages(t)
+        # State derivatives
+        rl_f = self.grid_filter.f(i_gs, u_cs, e_gs)
+        # List of state derivatives
+        return rl_f
+
+    def save(self, sol):
+        """
+        Save the solution.
+
+        Parameters
+        ----------
+        sol : Bunch object
+            Solution from the solver.
+
+        """
+        self.data.t.extend(sol.t)
+        self.data.i_gs.extend(sol.y[0])
+        self.data.q.extend(sol.q)
+
+
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format and post-process them.
+
+        """
+        # From lists to the ndarray
+        self.data.t = np.asarray(self.data.t)
+        self.data.i_gs = np.asarray(self.data.i_gs)
+        self.data.q = np.asarray(self.data.q)
+
+        # Some useful variables
+        self.data.e_gs = self.grid_model.voltages(self.data.t)
+        self.data.theta = np.mod(self.data.t*self.grid_model.w_N, 2*np.pi)
+        self.data.u_cs = self.converter.ac_voltage(self.data.q, self.converter.u_dc0)
+        self.data.u_gs = self.grid_filter.pcc_voltages(
+            self.data.i_gs,
+            self.data.u_cs,
+            self.data.e_gs)
+
+
+# %%
+class StiffSourceAndLCLFilterModel:
+    """
+    Continuous-time model for a stiff grid model with an LCL impedance model.
+
+    Parameters
+    ----------
+    grid_filter : LCLFilter
+        LCL dynamic model.
+    grid_model : StiffSource
+        Constant voltage source model.
+    converter : Inverter | PWMInverter
+        Inverter model.
+
+    """
+
+    def __init__(self, grid_filter=None, grid_model=None, converter=None):
+        self.grid_filter = grid_filter
+        self.grid_model = grid_model
+        self.converter = converter
+
+        # Initial time
+        self.t0 = 0
+
+        # Store the solution in these lists
+        self.data = Bunch()
+        self.data.t, self.data.q = [], []
+        self.data.i_gs = []
+        self.data.i_cs = []
+        self.data.u_fs = []
+
+    def get_initial_values(self):
+        """
+        Get the initial values.
+
+        Returns
+        -------
+        x0 : complex list, length 3
+            Initial values of the state variables.
+
+        """
+        x0 = [
+            self.grid_filter.i_cs0,
+            self.grid_filter.u_fs0,
+            self.grid_filter.i_gs0]
+
+        return x0
+
+    def set_initial_values(self, t0, x0):
+        """
+        Set the initial values.
+
+        Parameters
+        ----------
+        x0 : complex ndarray
+            Initial values of the state variables.
+
+        """
+        self.t0 = t0
+        self.grid_filter.i_cs0 = x0[0]
+        self.grid_filter.u_fs0 = x0[1]
+        self.grid_filter.i_gs0 = x0[2]
+        # calculation of grid-side voltage
+        e_gs0 = self.grid_model.voltages(t0)
+        # update pcc voltage
+        self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(
+                                                x0[2], x0[1], e_gs0)
+
+    def f(self, t, x):
+        """
+        Compute the complete state derivative list for the solver.
+
+        Parameters
+        ----------
+        t : float
+            Time.
+        x : complex ndarray
+            State vector.
+
+        Returns
+        -------
+        complex list
+            State derivatives.
+
+        """
+        # Unpack the states
+        i_cs, u_fs, i_gs = x
+        # Interconnections: outputs for computing the state derivatives
+        u_cs = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
+        e_gs = self.grid_model.voltages(t)
+        # State derivatives
+        lcl_f = self.grid_filter.f(i_cs, u_fs, i_gs, u_cs, e_gs)
+        # List of state derivatives
+        return lcl_f
+
+    def save(self, sol):
+        """
+        Save the solution.
+
+        Parameters
+        ----------
+        sol : Bunch object
+            Solution from the solver.
+
+        """
+        self.data.t.extend(sol.t)
+        self.data.i_cs.extend(sol.y[0])
+        self.data.u_fs.extend(sol.y[1])
+        self.data.i_gs.extend(sol.y[2])
+
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format and post-process them.
+
+        """
+        # From lists to the ndarray
+        self.data.t = np.asarray(self.data.t)
+        self.data.i_cs = np.asarray(self.data.i_cs)
+        self.data.u_fs = np.asarray(self.data.u_fs)
+        self.data.i_gs = np.asarray(self.data.i_gs)
+        self.data.u_dc = self.converter.u_dc0
+        self.data.q = np.asarray(self.data.q)
+        #self.data.theta = np.asarray(self.data.theta)
+        # Some useful variables
+        self.data.e_gs = self.grid_model.voltages(self.data.t)
+        self.data.theta = np.mod(self.data.t*self.grid_model.w_N, 2*np.pi)
+        self.data.u_cs = self.converter.ac_voltage(self.data.q, self.converter.u_dc0)
+        self.data.u_gs = self.grid_filter.pcc_voltages(
+            self.data.i_gs,
+            self.data.u_fs,
+            self.data.e_gs)
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
new file mode 100644
index 000000000..5b6c15210
--- /dev/null
+++ b/motulator/grid/model/_grid_filter.py
@@ -0,0 +1,298 @@
+"""
+Grid and AC filter impedance models.
+
+This module contains continuous-time models for subsystems comprising an AC 
+filter and a grid impedance between the converter and grid voltage sources. The 
+models are implemented with space vectors in stationary coordinates.
+
+"""
+import numpy as np
+from motulator.common.utils._utils import complex2abc
+
+
+# %%
+class LFilter:
+    """
+    Dynamic model for an inductive L filter and an inductive-resistive grid.
+
+    An L filter and an inductive-resistive grid impedance, between the converter 
+    and grid voltage sources, are modeled combining their inductances and series
+    resistances in a state equation. The grid current is used as a state 
+    variable. The point-of-common-coupling (PCC) voltage between the L filter 
+    and the grid impedance is separately calculated.
+
+    Parameters
+    ----------
+    L_f : float
+        Filter inductance (H)
+    R_f : float
+        Filter series resistance (Ω)
+    L_g : float
+        Grid inductance (H)
+    R_g : float
+        Grid resistance (Ω)
+
+    """
+    def __init__(self, U_gN=400*np.sqrt(2/3), L_f = 6e-3, R_f=0, L_g=0, R_g=0):
+        self.L_f = L_f
+        self.R_f = R_f
+        self.L_g = L_g
+        self.R_g = R_g
+        # Storing the PCC voltage value
+        self.u_gs0 = U_gN + 0j
+        # Initial values
+        self.i_gs0 = 0j
+
+
+    def pcc_voltages(self, i_gs, u_cs, e_gs):
+        """
+        Compute the PCC voltage between the L filter and grid impedance.
+        
+        Parameters
+        ----------
+        i_gs : complex
+            Grid current (A).
+        u_cs : complex
+            Converter voltage (V).
+        e_gs : complex
+            Grid voltage (V).
+
+        Returns
+        -------
+        u_gs : complex
+            Voltage at the point of common coupling (V).
+
+        """
+
+        # PCC voltage in stationary coordinates
+        u_gs = (self.L_g*u_cs + self.L_f*e_gs + 
+            (self.R_g*self.L_f - self.R_f*self.L_g)*i_gs)/(self.L_g+self.L_f)
+
+        return u_gs
+
+
+    def f(self, i_gs, u_cs, e_gs):
+        # pylint: disable=R0913
+        """
+        Compute the state derivatives.
+
+        Parameters
+        ----------
+        i_gs : complex
+            Grid current (A).
+        u_cs : complex
+            Converter voltage (V).
+        e_gs : complex
+            Grid voltage (V).
+
+        Returns
+        -------
+        complex list, length 1
+            Time derivative of the complex state vector, [di_gs]
+
+        """
+        # Calculation of the total impedance
+        L_t = self.L_f + self.L_g
+        R_t = self.R_f + self.R_g
+        di_gs = (u_cs - e_gs - R_t*i_gs)/L_t
+        return di_gs
+
+    def meas_currents(self):
+        """
+        Measure the phase currents at the end of the sampling period.
+
+        Returns
+        -------
+        i_g_abc : 3-tuple of floats
+            Grid phase currents (A).
+
+        """
+        # Grid phase currents from the corresponding space vector
+        i_g_abc = complex2abc(self.i_gs0)
+        return i_g_abc
+
+
+    def meas_pcc_voltage(self):
+        """
+        Measure the phase voltages at PCC at the end of the sampling period.
+
+        Returns
+        -------
+        u_g_abc : 3-tuple of floats
+            Phase voltage at the point of common coupling (V).
+
+        """  
+        # PCC phase voltages from the corresponding space vector
+        u_g_abc = complex2abc(self.u_gs0)
+        return u_g_abc
+
+
+    # %%
+class LCLFilter:
+    """
+    Dynamic model for an inductive-capacitive-inductive (LCL) filter and a grid.
+
+    An LCL filter and an inductive-resistive grid impedance, between the 
+    converter and grid voltage sources, are modeled using converter current, 
+    LCL-filter capacitor voltage and grid current as state variables. The grid 
+    inductance and resistance are included in the state equation of the grid 
+    current. The point-of-common-coupling (PCC) voltage between the LCL filter 
+    and the grid impedance is separately calculated.
+
+    Parameters
+    ----------
+    L_fc : float
+        Converter-side inductance of the LCL filter (H)
+    R_fc : float
+        Converter-side series resistance (Ω)
+    L_fg : float
+        Grid-side inductance of the LCL filter (H)
+    R_fg : float
+        Grid-side series resistance (Ω)
+    C_f : float
+        Capacitance of the LCL Filter (F)
+    G_f : float
+        Conductance of a resistor in parallel with the LCL-filter capacitor (S)
+    L_g : float
+        Grid inductance (H)
+    R_g : float
+        Grid resistance (Ω)
+
+
+    """
+    def __init__(
+            self, U_gN=400*np.sqrt(2/3), L_fc=6e-3, R_fc=0, L_fg=3e-3,
+            R_fg=0, C_f=10e-6, G_f=0, L_g=0, R_g=0):
+        self.L_fc = L_fc
+        self.R_fc = R_fc
+        self.L_fg = L_fg
+        self.R_fg = R_fg
+        self.C_f = C_f
+        self.G_f = G_f
+        self.L_g = L_g
+        self.R_g = R_g
+        # Storing useful variables
+        self.u_gs0 = U_gN + 0j
+        # Initial values
+        self.i_cs0 = 0j
+        self.i_gs0 = 0j
+        self.u_fs0 = U_gN + 0j
+
+
+    def pcc_voltages(self, i_gs, u_fs, e_gs):
+        """
+        Compute the PCC voltage between the LCL filter and the grid impedance.
+
+        Parameters
+        ----------
+        i_gs : complex
+            Grid current (A).
+        u_fs : complex
+            LCL-filter capacitor voltage (V).
+        e_gs : complex
+            Grid voltage (V).
+
+        Returns
+        -------
+        u_gs : complex
+            Voltage at the point of common coupling (V).
+
+        """
+        # PCC voltage in alpha-beta coordinates
+        u_gs = (self.L_fg*e_gs + self.L_g*u_fs + 
+            (self.R_g*self.L_fg-self.R_fg*self.L_g)*i_gs)/(self.L_g+self.L_fg)
+
+        return u_gs
+
+    def f(self, i_cs, u_fs, i_gs, u_cs, e_gs):
+        # pylint: disable=R0913
+        """
+        Compute the state derivatives.
+
+        Parameters
+        ----------
+        i_cs : complex
+            Converter current (A).
+        u_fs : complex
+            LCL-filter capacitor voltage (V).
+        i_gs : complex
+            Grid current (A).
+        u_cs : complex
+            Converter voltage (V).
+        e_gs : complex
+            Grid voltage (V).
+
+        Returns
+        -------
+        complex list, length 3
+            Time derivative of the complex state vector, [di_cs, du_fs, di_gs]
+
+        """
+        # Converter current dynamics
+        di_cs = (u_cs - u_fs - self.R_fc*i_cs)/self.L_fc
+        # Capacitor voltage dynamics
+        du_fs = (i_cs - i_gs - self.G_f*u_fs)/self.C_f
+        # Calculation of the total grid-side impedance
+        L_t = self.L_fg + self.L_g
+        R_t = self.R_fg + self.R_g
+        # Grid current dynamics
+        di_gs = (u_fs - e_gs - R_t*i_gs)/L_t
+
+        return [di_cs, du_fs, di_gs]
+
+    def meas_currents(self):
+        """
+        Measure the converter phase currents at the end of the sampling period.
+
+        Returns
+        -------
+        i_c_abc : 3-tuple of floats
+            Converter phase currents (A).
+
+        """
+        # Converter phase currents from the corresponding space vector
+        i_c_abc = complex2abc(self.i_cs0)
+
+        return i_c_abc
+
+    def meas_grid_currents(self):
+        """
+        Measure the grid phase currents at the end of the sampling period.
+
+        Returns
+        -------
+        i_g_abc : 3-tuple of floats
+            Grid phase currents (A).
+
+        """
+        # Grid phase currents from the corresponding space vector
+        i_g_abc = complex2abc(self.i_gs0)
+        return i_g_abc
+
+    def meas_cap_voltage(self):
+        """
+        Measure the capacitor phase voltages at the end of the sampling period.
+
+        Returns
+        -------
+        u_f_abc : 3-tuple of floats
+            Phase voltages of the LCL filter capacitor (V).
+
+        """  
+        # Capacitor phase voltages from the corresponding space vector
+        u_f_abc = complex2abc(self.u_fs0)
+        return u_f_abc
+
+    def meas_pcc_voltage(self):
+        """
+        Measure the PCC voltages at the end of the sampling period.
+
+        Returns
+        -------
+        u_g_abc : 3-tuple of floats
+            Phase voltages at the point of common coupling (V).
+
+        """
+        # PCC phase voltages from the corresponding space vector
+        u_g_abc = complex2abc(self.u_gs0)
+        return u_g_abc
diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
new file mode 100644
index 000000000..0f28764e2
--- /dev/null
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -0,0 +1,251 @@
+"""
+3-phase AC voltage source models.
+
+Two voltage sources with variable voltage magnitude, to simulate voltage dips
+and symmetrical short circuits are modeled. A stiff model with a constant
+frequency and a dynamic model with the electromechanical dynamics of a
+synchronous generator are considered. In this module, all space vectors are in
+stationary coordinates.
+
+"""
+
+import numpy as np
+from motulator.common.utils._utils import (complex2abc, abc2complex)
+
+# %%
+class StiffSource:
+    """
+    Grid subsystem.
+
+    This model is a constant frequency 3-phase voltage source of the AC grid.
+
+    Parameters
+    ----------
+    w_N : float
+        grid constant frequency (rad/s).
+    e_g_abs : function
+        3-phase grid voltage magnitude, phase-to-ground peak value (V).
+    """
+
+    def __init__(self, w_N=2*np.pi*50,
+                 e_g_abs=lambda t: 400*np.sqrt(2/3)):
+        self.w_N = w_N
+        self.e_g_abs = e_g_abs
+
+    def voltages(self, t):
+        """
+        Compute the grid voltage in stationary frame.
+           
+        Parameters
+        ----------
+        t : float
+            Time (s).
+
+        Returns
+        -------
+        e_gs: complex
+            grid complex voltage (V).
+
+        """
+        # Calculation of theta angle based on time and grid fixed frequency
+        theta = self.w_N*t
+
+        # Calculation of the three-phase voltage
+        e_g_a = self.e_g_abs(t)*np.cos(theta)
+        e_g_b = self.e_g_abs(t)*np.cos(theta-2*np.pi/3)
+        e_g_c = self.e_g_abs(t)*np.cos(theta-4*np.pi/3)
+
+
+        e_gs = abc2complex([e_g_a, e_g_b, e_g_c])
+        return e_gs
+
+
+    def meas_voltages(self, t):
+        """
+        Measure the phase voltages at the end of the sampling period.
+        
+        Parameters
+        ----------
+        t : float
+            Time (s).
+
+        Returns
+        -------
+        e_g_abc : 3-tuple of floats
+            Phase voltages (V).
+
+        """
+        # Grid voltage
+        e_g_abc = complex2abc(self.voltages(t))
+        return e_g_abc
+
+# %%
+class FlexSource:
+    """
+    Grid subsystem.
+    
+    This models the 3-phase voltage source of the AC grid while taking into
+    account the electromechanical dynamics of a typical grid generated by the 
+    synchronous generators.
+    
+    More information about the model can be found in [#ENT2013].
+    
+    [#ENT2013] : ENTSO-E, Documentation on Controller Tests in Test Grid
+    Configurations, Technical Report, 26.11.2013.
+    
+    Parameters
+    ----------
+    T_D : float
+        turbine delay time constant (s).
+    T_N : float
+        turbine derivative time constant (s).
+    H_g : float
+        grid inertia constant (s).
+    r_d : float
+        primary frequency droop control gain (p.u.).
+    T_gov : float
+        governor time constant (s).
+    w_N : float
+        grid constant frequency (rad/s).
+    S_grid : float
+        grid rated power (VA).
+    e_g_abs : function
+        3-phase grid voltage magnitude, phase-to-ground peak value (V).
+    p_m_ref : function
+        mechanical power output reference (W).
+    p_e : function
+        electrical power disturbance (W).
+        
+    """
+    def __init__(
+            self, T_D=10,
+            T_N=3,
+            H_g=3,
+            D_g=0,
+            r_d=.05,
+            T_gov=0.5,
+            w_N=2*np.pi*50,
+            S_grid =500e6,
+            e_g_abs=lambda t: 400*np.sqrt(2/3),
+            p_m_ref=lambda t: 0,
+            p_e=lambda t: 0):
+        self.e_g_abs=e_g_abs
+        self.w_N=w_N
+        self.S_grid=S_grid
+        self.T_D=T_D
+        self.T_N=T_N
+        self.H_g=H_g
+        self.D_g=D_g
+        self.r_d=r_d*w_N/S_grid
+        self.T_gov=T_gov
+        self.p_e=p_e
+        self.p_m_ref=p_m_ref
+        # Initial values
+        self.w_g0 = w_N
+        self.err_w_g0, self.p_gov0, self.x_turb0, self.theta_g0 = 0, 0, 0, 0
+
+    def f(self, t, err_w_g, p_gov, x_turb):
+        """
+        Compute the state derivatives.
+        
+        Parameters
+        ----------
+        t : float
+            Time (s).
+        err_w_g : float
+            grid angular speed deviation (mechanical rad/s).
+        p_gov : float
+            governor output power (W).
+        x_turb : float
+            turbine state variable (W).
+        p_e : float
+            electrical power disturbance (W).
+        Returns
+        -------
+        list, length 4
+            Time derivatives of the state vector.
+            
+        """
+        # calculation of mechanical power from the turbine output
+        p_m = (self.T_N/self.T_D)*p_gov + (1-(self.T_N/self.T_D))*x_turb
+        # swing equation
+        p_diff = (p_m - self.p_e(t))/self.S_grid # in per units
+        derr_w_g = self.w_N*(p_diff - self.D_g*err_w_g)/(2*self.H_g)
+        # governor dynamics
+        dp_gov = (self.p_m_ref(t) - (1/self.r_d)*err_w_g - p_gov) / self.T_gov
+        # turbine dynamics (lead-lag)
+        dx_turb = (p_gov - x_turb)/self.T_D
+        # integration of the angle
+        dtheta_g = self.w_N + err_w_g
+        return [derr_w_g, dp_gov, dx_turb, dtheta_g]
+
+    def voltages(self, t, theta_g):
+        """
+        Compute the grid voltage in stationary frame:
+           
+        Parameters
+        ----------
+        t : float
+            Time.
+        theta_g : float
+            grid electrical angle (rad).
+
+        Returns
+        -------
+        e_gs: complex
+            grid complex voltage (V).
+
+        """
+        # Calculation of the three-phase voltage
+        e_g_a = self.e_g_abs(t)*np.cos(theta_g)
+        e_g_b = self.e_g_abs(t)*np.cos(theta_g-2*np.pi/3)
+        e_g_c = self.e_g_abs(t)*np.cos(theta_g-4*np.pi/3)
+
+
+        e_gs = abc2complex([e_g_a, e_g_b, e_g_c])
+        return e_gs
+
+
+    def meas_voltages(self, t):
+        """
+        Measure the phase voltages at the end of the sampling period.
+
+        Returns
+        -------
+        e_g_abc : 3-tuple of floats
+            Phase voltages (V).
+
+        """
+        # Grid voltage
+        e_g_abc = complex2abc(self.voltages(t))
+        return e_g_abc
+
+
+    def meas_freq(self):
+        """
+        Measure the grid frequency.
+        
+        This returns the grid frequency at the end of the sampling period.
+        
+        Returns
+        -------
+        w_g0 : float
+            Grid angular speed (rad/s).
+            
+        """
+        w_g0 = self.w_g0
+        return w_g0
+
+    def meas_angle(self):
+        """
+        Measure the grid angle.
+        
+        This returns the grid angle at the end of the sampling period.
+        
+        Returns
+        -------
+        theta_g0 : float
+            grid electrical angle (rad).
+            
+        """
+        return self.theta_g0

From d99e681fe6681966a2172426323b865fdf114e28 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 09:53:25 +0300
Subject: [PATCH 05/26] added grid/control/_common.py

---
 .../grid/plot_grid_following_control_grid_converter_10kVA.py   | 1 -
 motulator/grid/control/_grid_following.py                      | 3 ++-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 7dd7de816..a3b934188 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -15,7 +15,6 @@
 import numpy as np
 from motulator.grid import model, control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
-add-gfl
 
 # To check the computation time of the program
 import time
diff --git a/motulator/grid/control/_grid_following.py b/motulator/grid/control/_grid_following.py
index 485e2718e..b0bd0ab0f 100644
--- a/motulator/grid/control/_grid_following.py
+++ b/motulator/grid/control/_grid_following.py
@@ -10,7 +10,8 @@
 import numpy as np
 
 from motulator.common.utils._utils import abc2complex
-from motulator.common.control import (Ctrl, PWM, ComplexFFPICtrl, Clock)
+from motulator.common.control import (PWM, ComplexFFPICtrl, Clock)
+from motulator.grid.control._common import (Ctrl, DCBusVoltCtrl)
 
 
 # %%

From 81afd2e6c53cbed019182975c6e08ba733f08215 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 09:56:17 +0300
Subject: [PATCH 06/26] added missing _common.py

---
 motulator/grid/control/_common.py | 90 +++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 motulator/grid/control/_common.py

diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
new file mode 100644
index 000000000..5d8394056
--- /dev/null
+++ b/motulator/grid/control/_common.py
@@ -0,0 +1,90 @@
+"""Common control functions and classes."""
+
+import numpy as np
+#from motulator.common.utils._utils import (abc2complex, complex2abc)
+from motulator.common.control._control import (Clock, PICtrl)
+from gritulator._utils import Bunch
+
+
+# %%
+class DCBusVoltCtrl(PICtrl):
+    """
+    PI DC-bus voltage controller.
+
+    This provides an interface for a DC-bus controller. The gains are
+    initialized based on the desired closed-loop bandwidth and the DC-bus
+    capacitance estimate. The PI controller is designed to control the energy
+    of the DC-bus capacitance and not the DC-bus voltage in order to have a
+    linear closed-loop system [#Hur2001]_.
+
+    Parameters
+    ----------
+    zeta : float
+        Damping ratio of the closed-loop system.
+    alpha_dc : float
+        Closed-loop bandwidth (rad/s). 
+    p_max : float, optional
+        Maximum converter power (W). The default is inf.
+        
+    References
+    ----------
+    .. [#Hur2001] Hur, Jung, Nam, "A Fast Dynamic DC-Link Power-Balancing
+       Scheme for a PWM Converter–Inverter System," IEEE Trans. Ind. Electron.,
+       2001, https://doi.org/10.1109/41.937412
+
+    """
+
+    def __init__(self, zeta, alpha_dc, p_max=np.inf):
+        k_p = -2*zeta*alpha_dc
+        k_i = -(alpha_dc**2)
+        k_t = k_p
+        super().__init__(k_p, k_i, k_t, p_max)
+
+
+# %%
+class Ctrl:
+    """Base class for the control system."""
+
+    def __init__(self):
+        self.data = Bunch()  # Data store
+        self.clock = Clock()  # Digital clock
+
+    def __call__(self, mdl):
+        """
+        Run the main control loop.
+
+        The main control loop is callable that returns the sampling
+        period `T_s` (float)  and the duty ratios `d_abc` (ndarray, shape (3,)) 
+        for the next sampling period.
+
+        Parameters
+        ----------
+        mdl : Model
+            System model containing methods for getting the feedback signals.
+
+        """
+        raise NotImplementedError
+
+    def save(self, data):
+        """
+        Save the internal date of the control system.
+
+        Parameters
+        ----------
+        data : bunch or dict
+            Contains the data to be saved.
+
+        """
+        for key, value in data.items():
+            self.data.setdefault(key, []).extend([value])
+
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format.
+
+        This method can be run after the simulation has been completed in order 
+        to simplify plotting and analysis of the stored data.
+
+        """
+        for key in self.data:
+            self.data[key] = np.asarray(self.data[key])

From 7cf7b6f8e416e2b3b09c450473c7306d9233cdb1 Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Mon, 10 Jun 2024 10:11:38 +0300
Subject: [PATCH 07/26] Changed Bunch datatypes to SimpleNamespace

---
 motulator/grid/control/_common.py         |  7 ++++---
 motulator/grid/model/_const_freq_model.py | 11 +++++------
 motulator/grid/utils/_plots.py            |  4 ++--
 3 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index 5d8394056..975a66346 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -3,7 +3,8 @@
 import numpy as np
 #from motulator.common.utils._utils import (abc2complex, complex2abc)
 from motulator.common.control._control import (Clock, PICtrl)
-from gritulator._utils import Bunch
+from types import SimpleNamespace
+#from gritulator._utils import Bunch
 
 
 # %%
@@ -46,7 +47,7 @@ class Ctrl:
     """Base class for the control system."""
 
     def __init__(self):
-        self.data = Bunch()  # Data store
+        self.data = SimpleNamespace()  # Data store
         self.clock = Clock()  # Digital clock
 
     def __call__(self, mdl):
@@ -71,7 +72,7 @@ def save(self, data):
 
         Parameters
         ----------
-        data : bunch or dict
+        data : SimpleNamespace or dict
             Contains the data to be saved.
 
         """
diff --git a/motulator/grid/model/_const_freq_model.py b/motulator/grid/model/_const_freq_model.py
index eca3d18c6..f24393905 100644
--- a/motulator/grid/model/_const_freq_model.py
+++ b/motulator/grid/model/_const_freq_model.py
@@ -10,8 +10,7 @@
 """
 
 import numpy as np
-
-#from gritulator._utils import Bunch
+from types import SimpleNamespace
 
 # %%
 class StiffSourceAndLFilterModel:
@@ -38,7 +37,7 @@ def __init__(self, grid_filter=None, grid_model=None, converter=None):
         self.t0 = 0
 
         # Store the solution in these lists
-        self.data = Bunch()
+        self.data = SimpleNamespace()
         self.data.t, self.data.q = [], []
         self.data.i_gs = []
 
@@ -109,7 +108,7 @@ def save(self, sol):
 
         Parameters
         ----------
-        sol : Bunch object
+        sol : SimpleNamespace object
             Solution from the solver.
 
         """
@@ -163,7 +162,7 @@ def __init__(self, grid_filter=None, grid_model=None, converter=None):
         self.t0 = 0
 
         # Store the solution in these lists
-        self.data = Bunch()
+        self.data = SimpleNamespace()
         self.data.t, self.data.q = [], []
         self.data.i_gs = []
         self.data.i_cs = []
@@ -239,7 +238,7 @@ def save(self, sol):
 
         Parameters
         ----------
-        sol : Bunch object
+        sol : SimpleNamespace object
             Solution from the solver.
 
         """
diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py
index b4a41b82b..75259cb08 100644
--- a/motulator/grid/utils/_plots.py
+++ b/motulator/grid/utils/_plots.py
@@ -6,7 +6,7 @@
 from cycler import cycler
 
 from motulator.common.utils import complex2abc
-from gritulator._utils import Bunch
+from types import SimpleNamespace
 
 # Plotting parameters
 plt.rcParams['axes.prop_cycle'] = cycler(color='brgcmyk')
@@ -57,7 +57,7 @@ def plot_grid(
     if base is None:
         pu_vals = False
         # Scaling with unity base values except for power use kW
-        base = Bunch(w=1, u=1, i=1, p=1000)
+        base = SimpleNamespace(w=1, u=1, i=1, p=1000)
     else:
         pu_vals = True
         

From 28d421a10d5c8e903a1d2fe90755b50f0f398d6e Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 11:09:29 +0300
Subject: [PATCH 08/26] changed imports for gfl example

---
 ..._following_control_grid_converter_10kVA.py | 11 ++++----
 motulator/grid/control/__init__.py            |  9 ++++++
 .../grid/control/grid_following/__init__.py   | 11 ++++++++
 .../{ => grid_following}/_grid_following.py   |  0
 motulator/grid/model/__init__.py              | 28 +++++++++++++++++++
 motulator/grid/model/_grid_filter.py          |  2 +-
 6 files changed, 55 insertions(+), 6 deletions(-)
 create mode 100644 motulator/grid/control/grid_following/__init__.py
 rename motulator/grid/control/{ => grid_following}/_grid_following.py (100%)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index a3b934188..7084eea7d 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -13,7 +13,8 @@
 # Imports.
 
 import numpy as np
-from motulator.grid import model, control
+from motulator.grid import model
+import motulator.grid.control.grid_following as control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
 
 # To check the computation time of the program
@@ -35,7 +36,7 @@
 # AC grid model (either constant frequency or dynamic electromechanical model)
 grid_model = model.StiffSource(w_N=2*np.pi*50)
 converter = model.Inverter(u_dc=650)
-mdl = model.ac_grid.StiffSourceAndLFilterModel(
+mdl = model.StiffSourceAndLFilterModel(
     grid_filter, grid_model, converter)
 
 
@@ -43,13 +44,13 @@
 # Configure the control system.
 
 # Control parameters
-pars = control.grid_following.GridFollowingCtrlPars(
+pars = control.GridFollowingCtrlPars(
             L_f=10e-3,
             f_sw = 5e3,
             T_s = 1/(10e3),
-            i_max = 1.5*base_values.i,
+            i_max = 1.5*base.i,
             )
-ctrl = control.grid_following.GridFollowingCtrl(pars)
+ctrl = control.GridFollowingCtrl(pars)
 
 
 # %%
diff --git a/motulator/grid/control/__init__.py b/motulator/grid/control/__init__.py
index 18bebfd8d..554d43115 100644
--- a/motulator/grid/control/__init__.py
+++ b/motulator/grid/control/__init__.py
@@ -1,2 +1,11 @@
 """Controls for grid-connected converters."""
 
+from motulator.grid.control._common import (
+    Ctrl,
+    DCBusVoltCtrl,
+)
+
+__all__ = [
+    "Ctrl",
+    "DCBusVoltCtrl"
+]
diff --git a/motulator/grid/control/grid_following/__init__.py b/motulator/grid/control/grid_following/__init__.py
new file mode 100644
index 000000000..bfa6cf012
--- /dev/null
+++ b/motulator/grid/control/grid_following/__init__.py
@@ -0,0 +1,11 @@
+"""This package contains example controllers for grid following converters."""
+
+from motulator.grid.control.grid_following._grid_following import (
+    GridFollowingCtrl,
+    GridFollowingCtrlPars,
+)
+
+__all__ = [
+    "GridFollowingCtrl",
+    "GridFollowingCtrlPars"
+]
diff --git a/motulator/grid/control/_grid_following.py b/motulator/grid/control/grid_following/_grid_following.py
similarity index 100%
rename from motulator/grid/control/_grid_following.py
rename to motulator/grid/control/grid_following/_grid_following.py
diff --git a/motulator/grid/model/__init__.py b/motulator/grid/model/__init__.py
index fd1e594ec..10b0ee2d8 100644
--- a/motulator/grid/model/__init__.py
+++ b/motulator/grid/model/__init__.py
@@ -5,7 +5,35 @@
     StiffSourceAndLCLFilterModel,
 )
 
+from motulator.grid.model._grid_filter import (
+    LFilter,
+    LCLFilter,
+)
+
+from motulator.grid.model._grid_volt_source import (
+    StiffSource,
+    FlexSource,
+)
+
+from motulator.common.model._converter import (
+    FrequencyConverter,
+    Inverter,
+)
+
+from motulator.common.model._simulation import (
+    CarrierComparison,
+    Simulation,
+)
+
 __all__ = [
     "StiffSourceAndLFilterModel",
     "StiffSourceAndLCLFilterModel",
+    "LFilter",
+    "LCLFilter",
+    "StiffSource",
+    "FlexSource",
+    "FrequencyConverter",
+    "Inverter",
+    "CarrierComparison",
+    "Simulation",
 ]
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index 5b6c15210..967f94645 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -127,7 +127,7 @@ def meas_pcc_voltage(self):
         return u_g_abc
 
 
-    # %%
+# %%
 class LCLFilter:
     """
     Dynamic model for an inductive-capacitive-inductive (LCL) filter and a grid.

From 71fe2bf150de4977be9e997362b35b1e417056b0 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 14:03:02 +0300
Subject: [PATCH 09/26] changed L-filter model to use Subsystem base class

---
 motulator/grid/model/_grid_filter.py | 82 +++++++++++-----------------
 1 file changed, 31 insertions(+), 51 deletions(-)

diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index 967f94645..7407d2621 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -6,12 +6,15 @@
 models are implemented with space vectors in stationary coordinates.
 
 """
+from types import SimpleNamespace
+
 import numpy as np
 from motulator.common.utils._utils import complex2abc
+from motulator.common.model import Subsystem
 
 
 # %%
-class LFilter:
+class LFilter(Subsystem):
     """
     Dynamic model for an inductive L filter and an inductive-resistive grid.
 
@@ -33,70 +36,47 @@ class LFilter:
         Grid resistance (Ω)
 
     """
-    def __init__(self, U_gN=400*np.sqrt(2/3), L_f = 6e-3, R_f=0, L_g=0, R_g=0):
-        self.L_f = L_f
-        self.R_f = R_f
-        self.L_g = L_g
-        self.R_g = R_g
-        # Storing the PCC voltage value
-        self.u_gs0 = U_gN + 0j
-        # Initial values
-        self.i_gs0 = 0j
-
+    def __init__(self, U_gN, L_f, R_f=0, L_g=0, R_g=0, u_gs=None):
+        super().__init__()
+        self.par = SimpleNamespace(U_gN=U_gN, L_f=L_f, R_f=R_f, L_g=L_g, R_g=R_g)
+        self._u_gs = u_gs
+        self.state = SimpleNamespace(i_gs=0)
+        self.sol_states = SimpleNamespace(i_gs=[])
 
-    def pcc_voltages(self, i_gs, u_cs, e_gs):
-        """
-        Compute the PCC voltage between the L filter and grid impedance.
-        
-        Parameters
-        ----------
-        i_gs : complex
-            Grid current (A).
-        u_cs : complex
-            Converter voltage (V).
-        e_gs : complex
-            Grid voltage (V).
-
-        Returns
-        -------
-        u_gs : complex
-            Voltage at the point of common coupling (V).
 
-        """
+    def u_gs(self):
+        """Compute the PCC voltage between the L filter and grid impedance."""
+        if callable(self._u_gs):
+            return self._u_gs(self.state.i_gs)
+        return (self.par.L_g*self.inp.u_cs + self.par.L_f*self.inp.e_gs +
+            (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
+            self.state.i_gs)/(self.par.L_g+self.par.L_f)
 
-        # PCC voltage in stationary coordinates
-        u_gs = (self.L_g*u_cs + self.L_f*e_gs + 
-            (self.R_g*self.L_f - self.R_f*self.L_g)*i_gs)/(self.L_g+self.L_f)
 
-        return u_gs
+    def set_outputs(self, _):
+        """Set output variables."""
+        state, out = self.state, self.out
+        out.i_gs = state.i_gs
 
 
-    def f(self, i_gs, u_cs, e_gs):
+    def rhs(self):
         # pylint: disable=R0913
         """
         Compute the state derivatives.
 
-        Parameters
-        ----------
-        i_gs : complex
-            Grid current (A).
-        u_cs : complex
-            Converter voltage (V).
-        e_gs : complex
-            Grid voltage (V).
-
         Returns
         -------
         complex list, length 1
             Time derivative of the complex state vector, [di_gs]
 
         """
-        # Calculation of the total impedance
-        L_t = self.L_f + self.L_g
-        R_t = self.R_f + self.R_g
-        di_gs = (u_cs - e_gs - R_t*i_gs)/L_t
+        state, inp, par = self.state, self.inp, self.par
+        L_t = par.L_f + par.L_g
+        R_t = par.R_f + par.R_g
+        di_gs = (inp.u_cs - inp.e_gs - R_t*state.i_gs)/L_t
         return di_gs
 
+
     def meas_currents(self):
         """
         Measure the phase currents at the end of the sampling period.
@@ -108,11 +88,11 @@ def meas_currents(self):
 
         """
         # Grid phase currents from the corresponding space vector
-        i_g_abc = complex2abc(self.i_gs0)
+        i_g_abc = complex2abc(self.state.i_gs)
         return i_g_abc
 
 
-    def meas_pcc_voltage(self):
+    def meas_voltages(self):
         """
         Measure the phase voltages at PCC at the end of the sampling period.
 
@@ -121,9 +101,9 @@ def meas_pcc_voltage(self):
         u_g_abc : 3-tuple of floats
             Phase voltage at the point of common coupling (V).
 
-        """  
+        """
         # PCC phase voltages from the corresponding space vector
-        u_g_abc = complex2abc(self.u_gs0)
+        u_g_abc = complex2abc(self.state.u_gs)
         return u_g_abc
 
 

From 932b9fd1cacfd41148c3303a71cbe543cc7959bd Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 14:14:35 +0300
Subject: [PATCH 10/26] added missing property attribute

---
 motulator/grid/model/_grid_filter.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index 7407d2621..8caf1caab 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -44,6 +44,7 @@ def __init__(self, U_gN, L_f, R_f=0, L_g=0, R_g=0, u_gs=None):
         self.sol_states = SimpleNamespace(i_gs=[])
 
 
+    @property
     def u_gs(self):
         """Compute the PCC voltage between the L filter and grid impedance."""
         if callable(self._u_gs):

From 9c486c797975de26204b505a98baa6ef34a8cf89 Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Mon, 10 Jun 2024 14:19:39 +0300
Subject: [PATCH 11/26] stiff source might work

---
 motulator/grid/control/_common.py         |  1 -
 motulator/grid/model/_grid_volt_source.py | 33 ++++++++++++-----------
 2 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index 975a66346..a4eb20085 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -4,7 +4,6 @@
 #from motulator.common.utils._utils import (abc2complex, complex2abc)
 from motulator.common.control._control import (Clock, PICtrl)
 from types import SimpleNamespace
-#from gritulator._utils import Bunch
 
 
 # %%
diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index 0f28764e2..249c4f1cc 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -10,10 +10,13 @@
 """
 
 import numpy as np
+
+from motulator.common.model import Subsystem
 from motulator.common.utils._utils import (complex2abc, abc2complex)
+from types import SimpleNamespace
 
 # %%
-class StiffSource:
+class StiffSource(Subsystem):
     """
     Grid subsystem.
 
@@ -29,8 +32,13 @@ class StiffSource:
 
     def __init__(self, w_N=2*np.pi*50,
                  e_g_abs=lambda t: 400*np.sqrt(2/3)):
-        self.w_N = w_N
-        self.e_g_abs = e_g_abs
+        super().__init__()
+        self.par = SimpleNamespace(w_N=w_N, e_g_abs=e_g_abs)
+        
+    @property
+    def w_N(self):
+        """Grid constant frequency (rad/s)."""
+        return self.par.w_N
 
     def voltages(self, t):
         """
@@ -80,7 +88,7 @@ def meas_voltages(self, t):
         return e_g_abc
 
 # %%
-class FlexSource:
+class FlexSource(Subsystem):
     """
     Grid subsystem.
     
@@ -129,17 +137,12 @@ def __init__(
             e_g_abs=lambda t: 400*np.sqrt(2/3),
             p_m_ref=lambda t: 0,
             p_e=lambda t: 0):
-        self.e_g_abs=e_g_abs
-        self.w_N=w_N
-        self.S_grid=S_grid
-        self.T_D=T_D
-        self.T_N=T_N
-        self.H_g=H_g
-        self.D_g=D_g
-        self.r_d=r_d*w_N/S_grid
-        self.T_gov=T_gov
-        self.p_e=p_e
-        self.p_m_ref=p_m_ref
+        
+        super().__init__()
+        self.par = SimpleNamespace(
+            T_D=T_D, T_N=T_N, H_g=H_g, D_g=D_g, r_d=r_d, T_gov=T_gov, w_N=w_N,
+            S_grid=S_grid, e_g_abs=e_g_abs, p_m_ref=p_m_ref, p_e=p_e)
+        
         # Initial values
         self.w_g0 = w_N
         self.err_w_g0, self.p_gov0, self.x_turb0, self.theta_g0 = 0, 0, 0, 0

From 06a0e3969a83dbca7c1666c524196448e56ee52b Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 14:59:41 +0300
Subject: [PATCH 12/26] Changed StiffSourceAndLFilterModel to use Model base
 class

---
 ..._following_control_grid_converter_10kVA.py |   5 +-
 motulator/grid/model/_const_freq_model.py     | 125 +++---------------
 motulator/grid/model/_grid_filter.py          |   6 +-
 3 files changed, 27 insertions(+), 109 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 7084eea7d..633c24b18 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -36,8 +36,7 @@
 # AC grid model (either constant frequency or dynamic electromechanical model)
 grid_model = model.StiffSource(w_N=2*np.pi*50)
 converter = model.Inverter(u_dc=650)
-mdl = model.StiffSourceAndLFilterModel(
-    grid_filter, grid_model, converter)
+mdl = model.StiffSourceAndLFilterModel(converter, grid_filter, grid_model)
 
 
 # %%
@@ -68,7 +67,7 @@
 # %%
 # Create the simulation object and simulate it.
 
-sim = model.Simulation(mdl, ctrl, pwm=False)
+sim = model.Simulation(mdl, ctrl)
 sim.simulate(t_stop = .1)
 
 # Print the execution time
diff --git a/motulator/grid/model/_const_freq_model.py b/motulator/grid/model/_const_freq_model.py
index f24393905..9cebd352f 100644
--- a/motulator/grid/model/_const_freq_model.py
+++ b/motulator/grid/model/_const_freq_model.py
@@ -8,12 +8,13 @@
     vector models are implemented in stationary coordinates.
 
 """
-
-import numpy as np
 from types import SimpleNamespace
+import numpy as np
+
+from motulator.common.model import Model
 
 # %%
-class StiffSourceAndLFilterModel:
+class StiffSourceAndLFilterModel(Model):
     """
     Continuous-time model for a stiff grid model with an RL impedance model.
 
@@ -28,113 +29,31 @@ class StiffSourceAndLFilterModel:
 
     """
 
-    def __init__(self, grid_filter=None, grid_model=None, converter=None):
+    def __init__(self, converter=None, grid_filter=None, grid_model=None):
+        super().__init__()
+        self.converter = converter
         self.grid_filter = grid_filter
         self.grid_model = grid_model
-        self.converter = converter
-
-        # Initial time
-        self.t0 = 0
-
-        # Store the solution in these lists
-        self.data = SimpleNamespace()
-        self.data.t, self.data.q = [], []
-        self.data.i_gs = []
+        self.subsystems = [self.converter, self.grid_filter, self.grid_model]
 
 
-    def get_initial_values(self):
-        """
-        Get the initial values.
-
-        Returns
-        -------
-        x0 : complex list, length 1
-            Initial values of the state variables.
-
-        """
-        x0 = [self.grid_filter.i_gs0]
-
-        return x0
-
-    def set_initial_values(self, t0, x0):
-        """
-        Set the initial values.
-
-        Parameters
-        ----------
-        x0 : complex ndarray
-            Initial values of the state variables.
-
-        """
-        self.t0 = t0
-        self.grid_filter.i_gs0 = x0[0]
-        # calculation of converter-side voltage
-        u_cs0 = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
-        # calculation of grid-side voltage
-        e_gs0 = self.grid_model.voltages(t0)
-        # update pcc voltage
-        self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(x0[0], u_cs0, e_gs0)
-
-    def f(self, t, x):
-        """
-        Compute the complete state derivative list for the solver.
-
-        Parameters
-        ----------
-        t : float
-            Time.
-        x : complex ndarray
-            State vector.
-
-        Returns
-        -------
-        complex list
-            State derivatives.
-
-        """
-        # Unpack the states
-        i_gs = x
-        # Interconnections: outputs for computing the state derivatives
-        u_cs = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
-        e_gs = self.grid_model.voltages(t)
-        # State derivatives
-        rl_f = self.grid_filter.f(i_gs, u_cs, e_gs)
-        # List of state derivatives
-        return rl_f
-
-    def save(self, sol):
-        """
-        Save the solution.
-
-        Parameters
-        ----------
-        sol : SimpleNamespace object
-            Solution from the solver.
-
-        """
-        self.data.t.extend(sol.t)
-        self.data.i_gs.extend(sol.y[0])
-        self.data.q.extend(sol.q)
+    def interconnect(self, _):
+        """Interconnect the subsystems."""
+        self.converter.inp.i_cs = self.grid_filter.out.i_cs
+        self.grid_filter.inp.u_cs = self.converter.out.u_cs
+        self.grid_filter.inp.e_gs = self.grid_model.out.e_gs
 
 
     def post_process(self):
-        """
-        Transform the lists to the ndarray format and post-process them.
-
-        """
-        # From lists to the ndarray
-        self.data.t = np.asarray(self.data.t)
-        self.data.i_gs = np.asarray(self.data.i_gs)
-        self.data.q = np.asarray(self.data.q)
-
-        # Some useful variables
-        self.data.e_gs = self.grid_model.voltages(self.data.t)
-        self.data.theta = np.mod(self.data.t*self.grid_model.w_N, 2*np.pi)
-        self.data.u_cs = self.converter.ac_voltage(self.data.q, self.converter.u_dc0)
-        self.data.u_gs = self.grid_filter.pcc_voltages(
-            self.data.i_gs,
-            self.data.u_cs,
-            self.data.e_gs)
+        """Post-process the solution."""
+        # Post-processing based on the states
+        super().post_process_states()
+        # Add the input data to the subsystems for post-processing
+        self.converter.data.i_cs = self.grid_filter.data.i_cs
+        self.grid_filter.data.u_cs = self.converter.data.u_cs
+        self.grid_filter.data.e_gs = self.grid_model.data.e_gs
+        # Post-processing based on the inputs and the states
+        super().post_process_with_inputs()
 
 
 # %%
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index 8caf1caab..354e42ae1 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -36,9 +36,9 @@ class LFilter(Subsystem):
         Grid resistance (Ω)
 
     """
-    def __init__(self, U_gN, L_f, R_f=0, L_g=0, R_g=0, u_gs=None):
+    def __init__(self, L_f, R_f=0, L_g=0, R_g=0, u_gs=None):
         super().__init__()
-        self.par = SimpleNamespace(U_gN=U_gN, L_f=L_f, R_f=R_f, L_g=L_g, R_g=R_g)
+        self.par = SimpleNamespace(L_f=L_f, R_f=R_f, L_g=L_g, R_g=R_g)
         self._u_gs = u_gs
         self.state = SimpleNamespace(i_gs=0)
         self.sol_states = SimpleNamespace(i_gs=[])
@@ -57,7 +57,7 @@ def u_gs(self):
     def set_outputs(self, _):
         """Set output variables."""
         state, out = self.state, self.out
-        out.i_gs = state.i_gs
+        out.i_gs, out.i_cs = state.i_gs, state.i_gs
 
 
     def rhs(self):

From 8d67072b231f6f1b5cbf927eaf88825d3cdbdf6c Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Mon, 10 Jun 2024 15:29:59 +0300
Subject: [PATCH 13/26] fixing stuff here and there

---
 ...plot_grid_following_control_grid_converter_10kVA.py | 10 ++++++----
 motulator/grid/control/__init__.py                     |  2 +-
 motulator/grid/utils/_helpers.py                       |  2 --
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 7084eea7d..188a66a45 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -13,8 +13,9 @@
 # Imports.
 
 import numpy as np
-from motulator.grid import model
-import motulator.grid.control.grid_following as control
+from motulator.grid import model, control
+import motulator.grid.control.grid_following as control   
+#import motulator.grid.control.grid_following as control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
 
 # To check the computation time of the program
@@ -32,12 +33,13 @@
 # Create the system model.
 
 # grid impedance and filter model
-grid_filter = model.LFilter(L_f=10e-3, L_g=0, R_g=0)
+grid_filter = model.LFilter(U_gN=400*np.sqrt(2/3) ,R_f=0 ,L_f=10e-3, L_g=0, R_g=0)
 # AC grid model (either constant frequency or dynamic electromechanical model)
 grid_model = model.StiffSource(w_N=2*np.pi*50)
 converter = model.Inverter(u_dc=650)
 mdl = model.StiffSourceAndLFilterModel(
     grid_filter, grid_model, converter)
+mdl.pwm = None  # disable the PWM model
 
 
 # %%
@@ -68,7 +70,7 @@
 # %%
 # Create the simulation object and simulate it.
 
-sim = model.Simulation(mdl, ctrl, pwm=False)
+sim = model.Simulation(mdl, ctrl)
 sim.simulate(t_stop = .1)
 
 # Print the execution time
diff --git a/motulator/grid/control/__init__.py b/motulator/grid/control/__init__.py
index 554d43115..69b108555 100644
--- a/motulator/grid/control/__init__.py
+++ b/motulator/grid/control/__init__.py
@@ -7,5 +7,5 @@
 
 __all__ = [
     "Ctrl",
-    "DCBusVoltCtrl"
+    "DCBusVoltCtrl",
 ]
diff --git a/motulator/grid/utils/_helpers.py b/motulator/grid/utils/_helpers.py
index f7db47a21..c7443e6a5 100644
--- a/motulator/grid/utils/_helpers.py
+++ b/motulator/grid/utils/_helpers.py
@@ -64,9 +64,7 @@ class BaseValues:
     i: float
     w: float
     psi: float
-    f: float
     p: float
-    f: float
     Z: float
     L: float
 

From 866b6b045ae813d996586b11c760ac2fee9f0958 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 10 Jun 2024 15:37:18 +0300
Subject: [PATCH 14/26] edit grid voltage source

---
 motulator/grid/model/_grid_volt_source.py | 18 ++++++++++++------
 motulator/grid/utils/_helpers.py          |  2 --
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index 249c4f1cc..da47a6a75 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -8,12 +8,12 @@
 stationary coordinates.
 
 """
+from types import SimpleNamespace
 
 import numpy as np
 
 from motulator.common.model import Subsystem
 from motulator.common.utils._utils import (complex2abc, abc2complex)
-from types import SimpleNamespace
 
 # %%
 class StiffSource(Subsystem):
@@ -34,12 +34,14 @@ def __init__(self, w_N=2*np.pi*50,
                  e_g_abs=lambda t: 400*np.sqrt(2/3)):
         super().__init__()
         self.par = SimpleNamespace(w_N=w_N, e_g_abs=e_g_abs)
-        
+
+
     @property
     def w_N(self):
         """Grid constant frequency (rad/s)."""
         return self.par.w_N
 
+
     def voltages(self, t):
         """
         Compute the grid voltage in stationary frame.
@@ -59,15 +61,19 @@ def voltages(self, t):
         theta = self.w_N*t
 
         # Calculation of the three-phase voltage
-        e_g_a = self.e_g_abs(t)*np.cos(theta)
-        e_g_b = self.e_g_abs(t)*np.cos(theta-2*np.pi/3)
-        e_g_c = self.e_g_abs(t)*np.cos(theta-4*np.pi/3)
-
+        e_g_a = self.par.e_g_abs(t)*np.cos(theta)
+        e_g_b = self.par.e_g_abs(t)*np.cos(theta-2*np.pi/3)
+        e_g_c = self.par.e_g_abs(t)*np.cos(theta-4*np.pi/3)
 
         e_gs = abc2complex([e_g_a, e_g_b, e_g_c])
         return e_gs
 
 
+    def set_outputs(self, t):
+        """Set output variables."""
+        self.out.e_gs = self.voltages(t)
+
+
     def meas_voltages(self, t):
         """
         Measure the phase voltages at the end of the sampling period.
diff --git a/motulator/grid/utils/_helpers.py b/motulator/grid/utils/_helpers.py
index f7db47a21..c7443e6a5 100644
--- a/motulator/grid/utils/_helpers.py
+++ b/motulator/grid/utils/_helpers.py
@@ -64,9 +64,7 @@ class BaseValues:
     i: float
     w: float
     psi: float
-    f: float
     p: float
-    f: float
     Z: float
     L: float
 

From e08f0da4c742feeda3de661fb134349a38aead0c Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Tue, 11 Jun 2024 13:23:48 +0300
Subject: [PATCH 15/26] initial adding of GridConverterCtrl base class

---
 motulator/grid/control/_common.py             | 146 ++++++++++++++----
 .../control/grid_following/_grid_following.py |   2 +-
 2 files changed, 116 insertions(+), 32 deletions(-)

diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index a4eb20085..bfd98584f 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -1,8 +1,9 @@
 """Common control functions and classes."""
-
+from abc import ABC
 import numpy as np
 #from motulator.common.utils._utils import (abc2complex, complex2abc)
-from motulator.common.control._control import (Clock, PICtrl)
+from motulator.common.control import (Ctrl, PICtrl)
+from motulator.common.utils import abc2complex, wrap
 from types import SimpleNamespace
 
 
@@ -42,49 +43,132 @@ def __init__(self, zeta, alpha_dc, p_max=np.inf):
 
 
 # %%
-class Ctrl:
-    """Base class for the control system."""
+class GridConverterCtrl(Ctrl, ABC):
+    """
+    Base class for control of grid-connected converters.
+    """
+    def __init__(self, par, T_s, sensorless):
+        super().__init__(T_s)
+        self.par = par
+        self.sensorless = sensorless
+        self.speed_ctrl = None
+        self.observer = None
+        self.ref = SimpleNamespace()
+
+    def get_electrical_measurements(self, fbk, mdl):
+        """
+        Measure the currents and voltages.
+        
+        Parameters
+        ----------
+        fbk : SimpleNamespace
+            Measured signals are added to this object.
+        mdl : Model
+            Continuous-time system model.
 
-    def __init__(self):
-        self.data = SimpleNamespace()  # Data store
-        self.clock = Clock()  # Digital clock
+        Returns
+        -------
+        fbk : SimpleNamespace
+            Measured signals, containing the following fields:
+
+                u_dc : float
+                    DC-bus voltage (V).
+                i_ss : complex
+                    Stator current (A) in stator coordinates.
+                u_ss : complex
+                    Realized stator voltage (V) in stator coordinates. This
+                    signal is obtained from the PWM.
 
-    def __call__(self, mdl):
         """
-        Run the main control loop.
+        fbk.u_dc = mdl.converter.meas_dc_voltage()
+        fbk.i_ss = abc2complex(mdl.machine.meas_currents())
+        fbk.u_ss = self.pwm.get_realized_voltage()
 
-        The main control loop is callable that returns the sampling
-        period `T_s` (float)  and the duty ratios `d_abc` (ndarray, shape (3,)) 
-        for the next sampling period.
+        return fbk
 
+    def get_mechanical_measurements(self, fbk, mdl):
+        """
+        Measure the speed and position.
+        
         Parameters
         ----------
+        fbk : SimpleNamespace
+            Measured signals are added to this object.
         mdl : Model
-            System model containing methods for getting the feedback signals.
-
+            Continuous-time system model.
+
+        Returns
+        -------
+        fbk : SimpleNamespace
+            Measured signals, containing the following fields:
+
+                w_m : float
+                    Rotor speed (electrical rad/s).
+                theta_m : float
+                    Rotor position (electrical rad).
+    
         """
-        raise NotImplementedError
+        fbk.w_m = self.par.n_p*mdl.mechanics.meas_speed()
+        fbk.theta_m = wrap(self.par.n_p*mdl.mechanics.meas_position())
 
-    def save(self, data):
-        """
-        Save the internal date of the control system.
+        return fbk
 
-        Parameters
-        ----------
-        data : SimpleNamespace or dict
-            Contains the data to be saved.
+    def get_feedback_signals(self, mdl):
+        """Get the feedback signals."""
+        fbk = super().get_feedback_signals(mdl)
+        fbk = self.get_electrical_measurements(fbk, mdl)
+        if not self.sensorless:
+            fbk = self.get_mechanical_measurements(fbk, mdl)
+        if self.observer:
+            fbk = self.observer.output(fbk)
 
-        """
-        for key, value in data.items():
-            self.data.setdefault(key, []).extend([value])
+        return fbk
 
-    def post_process(self):
+    def get_torque_reference(self, fbk, ref):
         """
-        Transform the lists to the ndarray format.
+        Get the torque reference in vector control.
+
+        This method can be used in vector control to get the torque reference 
+        from the speed controller. If the speed controller method `speed_ctrl` 
+        is None, the torque reference is obtained directly from the reference.
 
-        This method can be run after the simulation has been completed in order 
-        to simplify plotting and analysis of the stored data.
+        Parameters
+        ----------
+        fbk : SimpleNamespace
+            Feedback signals. In speed-control mode, the measured or estimated
+            rotor speed `w_m` is used to compute the torque reference.
+        ref : SimpleNamespace
+            Reference signals, containing the digital time `t`. The speed and 
+            torque references are added to this object.
+
+        Returns
+        -------
+        ref : SimpleNamespace
+            Reference signals, containing the following fields:
+                
+                w_m : float
+                    Speed reference (electrical rad/s).
+                tau_M : float
+                    Torque reference (Nm).  
 
         """
-        for key in self.data:
-            self.data[key] = np.asarray(self.data[key])
+        if self.speed_ctrl:
+            # Speed-control mode
+            ref.w_m = self.ref.w_m(ref.t)
+            ref_w_M = ref.w_m/self.par.n_p
+            w_M = fbk.w_m/self.par.n_p
+            ref.tau_M = self.speed_ctrl.output(ref_w_M, w_M)
+        else:
+            # Torque-control mode
+            ref.w_m = None
+            ref.tau_M = self.ref.tau_M(ref.t)
+
+        return ref
+
+    def update(self, fbk, ref):
+        """Extend the base class method."""
+        super().update(fbk, ref)
+        if self.speed_ctrl:
+            self.speed_ctrl.update(ref.T_s, ref.tau_M)
+        if self.observer:
+            self.observer.update(ref.T_s, fbk)
diff --git a/motulator/grid/control/grid_following/_grid_following.py b/motulator/grid/control/grid_following/_grid_following.py
index b0bd0ab0f..60a851d75 100644
--- a/motulator/grid/control/grid_following/_grid_following.py
+++ b/motulator/grid/control/grid_following/_grid_following.py
@@ -67,7 +67,7 @@ class GridFollowingCtrlPars:
 
 
 # %%
-class GridFollowingCtrl(Ctrl):
+class GridFollowingCtrl(GridConverterCtrl):
     """
     Grid following control for power converters.
 

From 794c12f09ce461ebff35b9211d743e421eb193d6 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Tue, 11 Jun 2024 17:41:53 +0300
Subject: [PATCH 16/26] changed to use GFL control from gritulator side, fixes
 to L-filter PCC voltage

---
 ..._following_control_grid_converter_10kVA.py |   9 +-
 motulator/grid/control/__init__.py            |  15 +-
 motulator/grid/control/_common.py             | 605 +++++++++++++++++-
 .../control/grid_following/_grid_following.py |   6 +-
 motulator/grid/model/_grid_filter.py          |  24 +-
 5 files changed, 632 insertions(+), 27 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 367787b03..a25ca2461 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -12,14 +12,15 @@
 # %%
 # Imports.
 
+import time
 import numpy as np
-from motulator.grid import model, control
-import motulator.grid.control.grid_following as control   
+
+from motulator.grid import model
+import motulator.grid.control.grid_following as control
 #import motulator.grid.control.grid_following as control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
 
 # To check the computation time of the program
-import time
 start_time = time.time()
 
 
@@ -39,8 +40,8 @@
 converter = model.Inverter(u_dc=650)
 
 
-mdl.pwm = None  # disable the PWM model
 mdl = model.StiffSourceAndLFilterModel(converter, grid_filter, grid_model)
+mdl.pwm = None  # disable the PWM model
 
 
 
diff --git a/motulator/grid/control/__init__.py b/motulator/grid/control/__init__.py
index 69b108555..b17c8ad41 100644
--- a/motulator/grid/control/__init__.py
+++ b/motulator/grid/control/__init__.py
@@ -1,11 +1,22 @@
 """Controls for grid-connected converters."""
 
 from motulator.grid.control._common import (
-    Ctrl,
+    ComplexFFPICtrl,
+    ComplexPICtrl,
+    RateLimiter,
     DCBusVoltCtrl,
+    PICtrl,
+    PWM
 )
 
+import motulator.grid.control.grid_following
+
 __all__ = [
-    "Ctrl",
+    "ComplexFFPICtrl",
+    "ComplexPICtrl",
+    "RateLimiter",
     "DCBusVoltCtrl",
+    "PICtrl",
+    "PWM",
+    "grid_following",
 ]
diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index bfd98584f..1573191cc 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -1,10 +1,280 @@
 """Common control functions and classes."""
 from abc import ABC
+from types import SimpleNamespace
+
 import numpy as np
 #from motulator.common.utils._utils import (abc2complex, complex2abc)
-from motulator.common.control import (Ctrl, PICtrl)
-from motulator.common.utils import abc2complex, wrap
-from types import SimpleNamespace
+#from motulator.common.control import (Ctrl, PICtrl)
+from motulator.common.utils import (complex2abc, abc2complex, wrap)
+
+# %%
+class PWM:
+    """
+    Duty ratios and realized voltage for three-phase PWM.
+
+    This contains the computation of the duty ratios and the realized voltage. 
+    The realized voltage is computed based on the measured DC-bus voltage and 
+    the duty ratios. The digital delay effects are taken into account in the 
+    realized voltage, assuming the delay of a single sampling period. The total
+    delay is 1.5 sampling periods due to the PWM (or ZOH) delay [#Bae2003]_.
+
+    Parameters
+    ----------
+    six_step: bool, optional
+        Enable six-step operation in overmodulation. The default is False.
+    
+    Attributes
+    ----------
+    realized_voltage : complex
+        Realized voltage (V) in synchronous coordinates.
+
+    References
+    ----------
+    .. [#Bae2003] Bae, Sul, "A compensation method for time delay of 
+       full-digital synchronous frame current regulator of PWM AC drives," IEEE
+       Trans. Ind. Appl., 2003, https://doi.org/10.1109/TIA.2003.810660
+    
+    """
+
+    def __init__(self, six_step=False):
+        self.six_step = six_step
+        self.realized_voltage = 0
+        self._u_ref_lim_old = 0
+
+    @staticmethod
+    def six_step_overmodulation(u_s_ref, u_dc):
+        """
+        Overmodulation up to six-step operation.
+
+        This method modifies the angle of the voltage reference vector in the
+        overmodulation region such that the six-step operation is reached 
+        [#Bol1997]_.
+
+        Parameters
+        ----------
+        u_s_ref : complex
+            Reference voltage (V) in stator coordinates.
+        u_dc : float
+            DC-bus voltage (V).
+
+        Returns
+        -------
+        u_s_ref_mod : complex
+            Reference voltage (V) in stator coordinates.
+
+        References
+        ----------
+        .. [#Bol1997] Bolognani, Zigliotto, "Novel digital continuous control of 
+           SVM inverters in the overmodulation range," IEEE Trans. Ind. Appl.,
+           1997, https://doi.org/10.1109/28.568019
+
+        """
+        # Limited magnitude
+        r = np.min([np.abs(u_s_ref), 2/3*u_dc])
+
+        if np.sqrt(3)*r > u_dc:
+            # Angle and sector of the reference vector
+            theta = np.angle(u_s_ref)
+            sector = np.floor(3*theta/np.pi)
+
+            # Angle reduced to the first sector (at which sector == 0)
+            theta0 = theta - sector*np.pi/3
+
+            # Intersection angle, see Eq. (9)
+            alpha_g = np.pi/6 - np.arccos(u_dc/(np.sqrt(3)*r))
+
+            # Modify the angle according to Eq. (4)
+            if alpha_g <= theta0 <= np.pi/6:
+                theta0 = alpha_g
+            elif np.pi/6 <= theta0 <= np.pi/3 - alpha_g:
+                theta0 = np.pi/3 - alpha_g
+
+            # Modified reference voltage
+            u_s_ref_mod = r*np.exp(1j*(theta0 + sector*np.pi/3))
+        else:
+            u_s_ref_mod = u_s_ref
+
+        return u_s_ref_mod
+
+    @staticmethod
+    def duty_ratios(u_s_ref, u_dc):
+        """
+        Compute the duty ratios for three-phase PWM.
+
+        This computes the duty ratios using a symmetrical suboscillation
+        method. This method is identical to the standard space-vector PWM.
+
+        Parameters
+        ----------
+        u_s_ref : complex
+            Voltage reference in stator coordinates (V).
+        u_dc : float
+            DC-bus voltage (V).
+
+        Returns
+        -------
+        d_abc : ndarray, shape (3,)
+            Duty ratios.
+
+        """
+        # Phase voltages without the zero-sequence voltage
+        u_abc = complex2abc(u_s_ref)
+
+        # Symmetrization by adding the zero-sequence voltage
+        u_0 = .5*(np.amax(u_abc) + np.amin(u_abc))
+        u_abc -= u_0
+
+        # Preventing overmodulation by means of a minimum phase error method
+        m = (2./u_dc)*np.amax(u_abc)
+        if m > 1:
+            u_abc = u_abc/m
+
+        # Duty ratios
+        d_abc = .5 + u_abc/u_dc
+
+        return d_abc
+
+    # pylint: disable=too-many-arguments
+    def __call__(self, T_s, u_ref, u_dc, theta, w):
+        """
+        Compute the duty ratios and update the state.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+        u_ref : complex
+            Voltage reference in synchronous coordinates (V).
+        u_dc : float
+            DC-bus voltage.
+        theta : float
+            Angle of synchronous coordinates (rad).
+        w : float
+            Angular speed of synchronous coordinates (rad/s).
+
+        Returns
+        -------
+        d_abc : ndarray, shape (3,)
+            Duty ratios for the next sampling period.
+
+        """
+        # Advance the angle due to the computational delay (T_s) and the ZOH
+        # (PWM) delay (0.5*T_s)
+        theta_comp = theta + 1.5*T_s*w
+
+        # Voltage reference in stator coordinates
+        u_s_ref = np.exp(1j*theta_comp)*u_ref
+
+        # Modify angle in the overmodulation region
+        if self.six_step:
+            u_s_ref = self.six_step_overmodulation(u_s_ref, u_dc)
+
+        # Duty ratios
+        d_abc = self.duty_ratios(u_s_ref, u_dc)
+
+        # Realizable voltage
+        u_s_ref_lim = abc2complex(d_abc)*u_dc
+        u_ref_lim = np.exp(-1j*theta_comp)*u_s_ref_lim
+
+        # Update the voltage estimate for the next sampling instant.
+        self.realized_voltage = .5*(self._u_ref_lim_old + u_ref_lim)
+        self._u_ref_lim_old = u_ref_lim
+
+        return d_abc
+
+
+# %%
+class PICtrl:
+    """
+    2DOF PI controller.
+
+    This implements a discrete-time 2DOF PI controller, whose continuous-time
+    counterpart is::
+
+        u = k_t*y_ref - k_p*y + (k_i/s)*(y_ref - y)
+
+    where `u` is the controller output, `y_ref` is the reference signal, `y` is
+    the feedback signal, and `1/s` refers to integration. The standard PI
+    controller is obtained by choosing ``k_t = k_p``. The integrator anti-windup 
+    is implemented based on the realized controller output.
+
+    Notes
+    -----
+    This contoller can be used, e.g., as a speed controller. In this case, `y` 
+    corresponds to the rotor angular speed `w_M` and `u` to the torque reference 
+    `tau_M_ref`.
+
+    Parameters
+    ----------
+    k_p : float
+        Proportional gain.
+    k_i : float
+        Integral gain.
+    k_t : float, optional
+        Reference-feedforward gain. The default is `k_p`.
+    u_max : float, optional
+        Maximum controller output. The default is inf.
+
+    Attributes
+    ----------
+    v : float
+        Input disturbance estimate.
+    u_i : float
+        Integral state.
+
+    """
+
+    def __init__(self, k_p, k_i, k_t=None, u_max=np.inf):
+        self.u_max = u_max
+        # Gains
+        self.k_p = k_p
+        self.k_t = k_t if k_t is not None else k_p
+        self.alpha_i = k_i/self.k_t  # Inverse of the integration time T_i
+        # States
+        self.v, self.u_i = 0, 0
+
+    def output(self, y_ref, y):
+        """
+        Compute the controller output.
+
+        Parameters
+        ----------
+        y_ref : float
+            Reference signal.
+        y : float
+            Feedback signal.
+
+        Returns
+        -------
+        u : float
+            Controller output.
+
+        """
+        # Estimate of a disturbance input
+        self.v = self.u_i - (self.k_p - self.k_t)*y
+
+        # Controller output
+        u = self.k_t*(y_ref - y) + self.v
+
+        # Limit the controller output
+        u = np.clip(u, -self.u_max, self.u_max)
+
+        return u
+
+    def update(self, T_s, u_lim):
+        """
+        Update the integral state.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+        u_lim : float
+            Realized (limited) controller output. If the actuator does not
+            saturate, ``u_lim = u``.
+
+        """
+        self.u_i += T_s*self.alpha_i*(u_lim - self.v)
 
 
 # %%
@@ -42,6 +312,335 @@ def __init__(self, zeta, alpha_dc, p_max=np.inf):
         super().__init__(k_p, k_i, k_t, p_max)
 
 
+# %%
+class ComplexPICtrl:
+    """
+    2DOF synchronous-frame complex-vector PI controller.
+
+    This implements a discrete-time 2DOF synchronous-frame complex-vector PI
+    controller, whose continuous-time counterpart is [#Bri2000]_::
+
+        u = k_t*i_ref - k_p*i + (k_i + 1j*w*k_t)/s*(i_ref - i)
+
+    where `u` is the controller output, `i_ref` is the reference signal, `i` is
+    the feedback signal, `w` is the angular speed of synchronous coordinates, 
+    and `1/s` refers to integration. This algorithm is compatible with both real 
+    and complex signals. The 1DOF version is obtained by setting ``k_t = k_p``. 
+    The integrator anti-windup is implemented based on the realized controller 
+    output.
+
+    Parameters
+    ----------
+    k_p : float
+        Proportional gain.
+    k_i : float
+        Integral gain.
+    k_t : float, optional
+        Reference-feedforward gain. The default is `k_p`.
+
+    Attributes
+    ----------
+    v : complex
+        Input disturbance estimate.
+    u_i : complex
+        Integral state.
+
+    Notes
+    -----
+    This contoller can be used, e.g., as a current controller. In this case, `i`
+    corresponds to the stator current and `u` to the stator voltage.
+    
+    References
+    ----------
+    .. [#Bri2000] Briz, Degner, Lorenz, "Analysis and design of current 
+       regulators using complex vectors," IEEE Trans. Ind. Appl., 2000,
+       https://doi.org/10.1109/28.845057
+
+    """
+
+    def __init__(self, k_p, k_i, k_t=None):
+        # Gains
+        self.k_p = k_p
+        self.k_t = k_t if k_t is not None else k_p
+        self.alpha_i = k_i/self.k_t  # Inverse of the integration time T_i
+        # States
+        self.v, self.u_i = 0, 0
+
+    def output(self, i_ref, i):
+        """
+        Compute the controller output.
+
+        Parameters
+        ----------
+        i_ref : complex
+            Reference signal.
+        i : complex
+            Feedback signal.
+
+        Returns
+        -------
+        u : complex
+            Controller output.
+
+        """
+        # Disturbance input estimate
+        self.v = self.u_i - (self.k_p - self.k_t)*i
+
+        # Controller output
+        u = self.k_t*(i_ref - i) + self.v
+
+        return u
+
+    def update(self, T_s, u_lim, w):
+        """
+        Update the integral state.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+        u_lim : complex
+            Realized (limited) controller output. If the actuator does not
+            saturate, ``u_lim = u``.
+        w : float
+            Angular speed of the reference frame (rad/s). 
+
+        """
+        self.u_i += T_s*(self.alpha_i + 1j*w)*(u_lim - self.v)
+
+
+# %%
+class ComplexFFPICtrl:
+    """
+    2DOF Synchronous-frame complex-vector PI controller with feedforward.
+
+    This implements a discrete-time 2DOF synchronous-frame complex-vector PI
+    controller similar to [#Bri2000]_, with an additional feedforward signal.                        
+    The gain selection corresponding to internal-model-control (IMC) is
+    equivalent to the continuous-time version given in [#Har2009]_::
+
+        u = k_p*(i_ref - i) + k_i/s*(i_ref - i) + 1j*w*L_f*i + u_g_ff 
+                
+    where `u` is the controller output, `i_ref` is the reference signal, `i` is
+    the feedback signal, u_g_ff is the (filtered) feedforward signal, `w` is
+    the angular speed of synchronous coordinates, '1j*w*L_f' is the decoupling
+    term estimate, and `1/s` refers to integration. This algorithm is
+    compatible with both real and complex signals. The integrator anti-windup
+    is implemented based on the realized controller output.
+
+    Parameters
+    ----------
+    k_p : float
+        Proportional gain.
+    k_i : float
+        Integral gain.
+    k_t : float, optional
+        Reference-feedforward gain. The default is `k_p`.
+    L_f : float, optional
+        Synchronous frame decoupling gain. The default is 0.
+
+    Attributes
+    ----------
+    v : complex
+        Input disturbance estimate.
+    u_i : complex
+        Integral state.
+
+    Notes
+    -----
+    This contoller can be used, e.g., as a current controller. In this case,
+    `i` corresponds to the converter current and `u` to the converter voltage.
+    
+    References
+    ----------
+       
+    .. [#Har2009] Harnefors, Bongiorno, "Current controller design
+       for passivity of the input admittance," 2009 13th European Conference
+       on Power Electronics and Applications, Barcelona, Spain, 2009.
+
+    """
+
+    def __init__(self, k_p, k_i, k_t=None, L_f=None):
+        # Gains
+        self.k_p = k_p
+        self.k_t = k_t if k_t is not None else k_p
+        self.alpha_i = k_i/self.k_t  # Inverse of the integration time T_i
+        self.L_f = L_f if L_f is not None else 0
+        # States
+        self.v, self.u_i = 0, 0
+
+    def output(self, i_ref, i, u_ff, w):
+        """
+        Compute the controller output.
+
+        Parameters
+        ----------
+        i_ref : complex
+            Reference signal.
+        i : complex
+            Feedback signal.
+        u_ff : complex
+            Feedforward signal.
+        w : float
+            Angular speed of the reference frame (rad/s). 
+
+        Returns
+        -------
+        u : complex
+            Controller output.
+
+        """
+        # Disturbance input estimate
+        self.v = self.u_i - (self.k_p - self.k_t)*i + (u_ff
+        + 1j*w*self.L_f*i)
+
+        # Controller output
+        u = self.k_t*(i_ref - i) + self.v
+
+        return u
+
+    def update(self, T_s, u_lim):
+        """
+        Update the integral state.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+        u_lim : complex
+            Realized (limited) controller output. If the actuator does not
+            saturate, ``u_lim = u``.
+
+        """
+        self.u_i += T_s*self.alpha_i*(u_lim - self.v)
+
+
+# %%
+class RateLimiter:
+    """
+    Rate limiter.
+
+    Parameters
+    ----------
+    rate_limit : float, optional
+        Rate limit. The default is inf.
+
+    """
+
+    def __init__(self, rate_limit=np.inf):
+        self.rate_limit = rate_limit
+        self._y_old = 0
+
+    def __call__(self, T_s, u):
+        """
+        Limit the input signal.
+
+        Parameters
+        ----------
+        u : float
+            Input signal.
+
+        Returns
+        -------
+        T_s : float
+            Sampling period (s).
+        y : float
+            Rate-limited output signal.
+
+        Notes
+        -----
+        In this implementation, the falling rate limit equals the (negative)
+        rising rate limit. If needed, these limits can be separated with minor
+        modifications in the class.
+
+        """
+        rate = (u - self._y_old)/T_s
+
+        if rate > self.rate_limit:
+            # Limit rising rate
+            y = self._y_old + T_s*self.rate_limit
+        elif rate < -self.rate_limit:
+            # Limit falling rate
+            y = self._y_old - T_s*self.rate_limit
+        else:
+            y = u
+
+        # Store the limited output
+        self._y_old = y
+
+        return y
+
+
+# %%
+class Clock:
+    """Digital clock."""
+
+    def __init__(self):
+        self.t = 0
+        self.t_reset = 1e9
+
+    def update(self, T_s):
+        """
+        Update the digital clock.
+
+        Parameters
+        ----------
+        T_s : float
+            Sampling period (s).
+
+        """
+        self.t += T_s if self.t < self.t_reset else 0
+
+
+# %%
+class Ctrl:
+    """Base class for the control system."""
+
+    def __init__(self):
+        self.data = SimpleNamespace()  # Data store
+        self.clock = Clock()  # Digital clock
+
+    def __call__(self, mdl):
+        """
+        Run the main control loop.
+
+        The main control loop is callable that returns the sampling
+        period `T_s` (float)  and the duty ratios `d_abc` (ndarray, shape (3,)) 
+        for the next sampling period.
+
+        Parameters
+        ----------
+        mdl : Model
+            System model containing methods for getting the feedback signals.
+
+        """
+        raise NotImplementedError
+
+    def save(self, data):
+        """
+        Save the internal date of the control system.
+
+        Parameters
+        ----------
+        data : bunch or dict
+            Contains the data to be saved.
+
+        """
+        for key, value in data.items():
+            self.data.setdefault(key, []).extend([value])
+
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format.
+
+        This method can be run after the simulation has been completed in order 
+        to simplify plotting and analysis of the stored data.
+
+        """
+        for key in self.data:
+            self.data[key] = np.asarray(self.data[key])
+
+
 # %%
 class GridConverterCtrl(Ctrl, ABC):
     """
diff --git a/motulator/grid/control/grid_following/_grid_following.py b/motulator/grid/control/grid_following/_grid_following.py
index 60a851d75..1e016502a 100644
--- a/motulator/grid/control/grid_following/_grid_following.py
+++ b/motulator/grid/control/grid_following/_grid_following.py
@@ -10,8 +10,8 @@
 import numpy as np
 
 from motulator.common.utils._utils import abc2complex
-from motulator.common.control import (PWM, ComplexFFPICtrl, Clock)
-from motulator.grid.control._common import (Ctrl, DCBusVoltCtrl)
+#from motulator.common.control import (PWM, ComplexFFPICtrl, Clock)
+from motulator.grid.control._common import (Ctrl, PWM, Clock, ComplexFFPICtrl, DCBusVoltCtrl)
 
 
 # %%
@@ -67,7 +67,7 @@ class GridFollowingCtrlPars:
 
 
 # %%
-class GridFollowingCtrl(GridConverterCtrl):
+class GridFollowingCtrl(Ctrl):
     """
     Grid following control for power converters.
 
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index 354e42ae1..e23435fe8 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -36,28 +36,22 @@ class LFilter(Subsystem):
         Grid resistance (Ω)
 
     """
-    def __init__(self, L_f, R_f=0, L_g=0, R_g=0, u_gs=None):
+    def __init__(self, U_gN, L_f, R_f=0, L_g=0, R_g=0):
         super().__init__()
         self.par = SimpleNamespace(L_f=L_f, R_f=R_f, L_g=L_g, R_g=R_g)
-        self._u_gs = u_gs
+        #self.inp = SimpleNamespace(u_cs=0+0j)
+        self.out = SimpleNamespace(u_gs=U_gN+0j)  # Needed for direct feedthrough
         self.state = SimpleNamespace(i_gs=0)
         self.sol_states = SimpleNamespace(i_gs=[])
 
 
-    @property
-    def u_gs(self):
-        """Compute the PCC voltage between the L filter and grid impedance."""
-        if callable(self._u_gs):
-            return self._u_gs(self.state.i_gs)
-        return (self.par.L_g*self.inp.u_cs + self.par.L_f*self.inp.e_gs +
-            (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
-            self.state.i_gs)/(self.par.L_g+self.par.L_f)
-
-
     def set_outputs(self, _):
         """Set output variables."""
         state, out = self.state, self.out
-        out.i_gs, out.i_cs = state.i_gs, state.i_gs
+        u_gs = (self.par.L_g*self.inp.u_cs + self.par.L_f*self.inp.e_gs +
+            (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
+            self.state.i_gs)/(self.par.L_g+self.par.L_f)
+        out.i_gs, out.i_cs, out.u_gs = state.i_gs, state.i_gs, u_gs
 
 
     def rhs(self):
@@ -93,7 +87,7 @@ def meas_currents(self):
         return i_g_abc
 
 
-    def meas_voltages(self):
+    def meas_pcc_voltage(self):
         """
         Measure the phase voltages at PCC at the end of the sampling period.
 
@@ -104,7 +98,7 @@ def meas_voltages(self):
 
         """
         # PCC phase voltages from the corresponding space vector
-        u_g_abc = complex2abc(self.state.u_gs)
+        u_g_abc = complex2abc(self.out.u_gs)
         return u_g_abc
 
 

From a9592784d96f842308dad0c458ae5d500125adae Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Wed, 12 Jun 2024 09:28:36 +0300
Subject: [PATCH 17/26] data structure changes in GFL control

---
 motulator/grid/control/_common.py             | 19 ++++++++++++++-----
 .../control/grid_following/_grid_following.py |  9 +--------
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index 1573191cc..3f3d32592 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -616,7 +616,7 @@ def __call__(self, mdl):
         """
         raise NotImplementedError
 
-    def save(self, data):
+    def save(self, **kwargs):
         """
         Save the internal date of the control system.
 
@@ -626,8 +626,14 @@ def save(self, data):
             Contains the data to be saved.
 
         """
-        for key, value in data.items():
-            self.data.setdefault(key, []).extend([value])
+        for name, arg in kwargs.items():
+            if not hasattr(self.data, name):
+                setattr(self.data, name, SimpleNamespace())
+            data = getattr(self.data, name)
+            for key, value in vars(arg).items():
+                if not hasattr(data, key):
+                    setattr(data, key, [])
+                getattr(data, key).append(value)
 
     def post_process(self):
         """
@@ -637,8 +643,11 @@ def post_process(self):
         to simplify plotting and analysis of the stored data.
 
         """
-        for key in self.data:
-            self.data[key] = np.asarray(self.data[key])
+        self.data = SimpleNamespace()
+        for name, values in vars(self.data).items():
+            setattr(self.data, name, SimpleNamespace())
+            for key, value in vars(values).items():
+                setattr(getattr(self.data, name), key, np.asarray(value))
 
 
 # %%
diff --git a/motulator/grid/control/grid_following/_grid_following.py b/motulator/grid/control/grid_following/_grid_following.py
index 1e016502a..4f9e60419 100644
--- a/motulator/grid/control/grid_following/_grid_following.py
+++ b/motulator/grid/control/grid_following/_grid_following.py
@@ -202,14 +202,7 @@ def __call__(self, mdl):
                                            self.theta_p, self.w_g)
 
         # Data logging
-        data = SimpleNamespace(
-            w_c = w_pll, theta_c = self.theta_p, u_c_ref = u_c_ref,
-            u_c = u_c, i_c = i_c, abs_u_g = abs_u_g,
-            d_abc_ref = d_abc_ref, i_c_ref = i_c_ref, u_dc = u_dc,
-            t = self.clock.t, p_g_ref = p_g_ref, u_dc_ref = u_dc_ref,
-            q_g_ref=q_g_ref, u_g = u_g,
-                     )
-        self.save(data)
+        self.save(fbk=mdl.fbk, ref=mdl.ref)
 
         # Update the states
         self.theta_p = theta_pll

From 84fe7b84023ff6da177ed76250a0d8a631d964ca Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Wed, 12 Jun 2024 09:39:44 +0300
Subject: [PATCH 18/26] copied dc_model folder from gritulator and modified
 init & dc_bus.py

---
 motulator/grid/model/dc_bus/__init__.py      |  14 +
 motulator/grid/model/dc_bus/_dc_bus.py       |  96 ++++++
 motulator/grid/model/dc_bus/_dc_dyn_model.py | 311 +++++++++++++++++++
 3 files changed, 421 insertions(+)
 create mode 100644 motulator/grid/model/dc_bus/__init__.py
 create mode 100644 motulator/grid/model/dc_bus/_dc_bus.py
 create mode 100644 motulator/grid/model/dc_bus/_dc_dyn_model.py

diff --git a/motulator/grid/model/dc_bus/__init__.py b/motulator/grid/model/dc_bus/__init__.py
new file mode 100644
index 000000000..2b6655419
--- /dev/null
+++ b/motulator/grid/model/dc_bus/__init__.py
@@ -0,0 +1,14 @@
+"""Continuous-time DC-bus models."""
+
+from motulator.grid.model.dc_bus._dc_bus import DCBus
+
+from motulator.grid.model.dc_bus._dc_dyn_model import (
+    DCBusAndLFilterModel,
+    DCBusAndLCLFilterModel,
+)
+
+__all__ = [
+    "DCBus",
+    "DCBusAndLFilterModel",
+    "DCBusAndLCLFilterModel",
+]
diff --git a/motulator/grid/model/dc_bus/_dc_bus.py b/motulator/grid/model/dc_bus/_dc_bus.py
new file mode 100644
index 000000000..5417d2d17
--- /dev/null
+++ b/motulator/grid/model/dc_bus/_dc_bus.py
@@ -0,0 +1,96 @@
+"""
+Dynamic model of a DC bus.
+
+A DC bus between an external current source or sink and a converter is modeled 
+considering an equivalent circuit comprising a capacitor and parallel resistor.
+"""
+
+from motulator.common.utils import complex2abc
+
+# %%
+class DCBus:
+    """
+    DC bus model
+
+    This model is used to compute the capacitive DC bus dynamics. Dynamics are 
+    modeled with an equivalent circuit comprising a capacitor and its parallel 
+    resistor that is parametrized using a conductance. The capacitor voltage is 
+    used as a state variable.
+
+    Parameters
+    ----------
+    C_dc : float
+        DC-bus capacitance (F)
+    G_dc : float
+        Parallel conductance of the capacitor (S)
+    i_ext : function
+        External DC current, seen as disturbance, `i_ext(t)`.
+
+    """
+    def __init__(self, C_dc=1e-3, G_dc=0, i_ext=lambda t: 0, u_dc0 = 650):
+        self.C_dc = C_dc
+        self.G_dc = G_dc
+        self.i_ext = i_ext
+        # Initial values
+        self.u_dc0 = u_dc0
+    
+    @staticmethod
+    def dc_current(i_c_abc, q):
+        """
+        Compute the converter DC current from the switching states and phase 
+        currents.
+    
+        Parameters
+        ----------
+        i_c_abc : ndarray, shape (3,)
+            Phase currents (A).
+        q : complex ndarray, shape (3,)
+            Switching state vectors corresponding to the switching instants.
+            For example, the switching state q[1] is applied at the interval
+            t_n_sw[1].
+    
+        Returns
+        -------
+        i_dc: float
+            Converter DC current (A)
+    
+        """
+        # Duty ratio back into three-phase ratios
+        q_abc = complex2abc(q)
+        # Dot product
+        i_dc = q_abc[0]*i_c_abc[0] + q_abc[1]*i_c_abc[1] + q_abc[2]*i_c_abc[2]
+        return i_dc
+    
+    def rhs(self, t, u_dc, i_dc):
+        """
+        Compute the state derivatives.
+
+        Parameters
+        ----------
+        t : float
+            Time (s)
+        u_dc: float
+            DC bus voltage (V)
+        i_dc : float
+            Converter DC current (A)
+        Returns
+        -------
+        double list, length 1
+                Time derivative of the complex state vector, [du_dc]
+
+        """
+        # State derivative
+        du_dc = (self.i_ext(t) - i_dc - self.G_dc*u_dc)/self.C_dc
+        return du_dc
+
+    def meas_dc_voltage(self):
+        """
+        Measure the DC bus voltage at the end of the sampling period.
+    
+        Returns
+        -------
+        u_dc: float
+            DC bus voltage (V)
+    
+        """
+        return self.u_dc0
diff --git a/motulator/grid/model/dc_bus/_dc_dyn_model.py b/motulator/grid/model/dc_bus/_dc_dyn_model.py
new file mode 100644
index 000000000..bea8ca52d
--- /dev/null
+++ b/motulator/grid/model/dc_bus/_dc_dyn_model.py
@@ -0,0 +1,311 @@
+"""
+DC-bus and filter model interconnectors.
+
+    This module interconnects the subsystems of a converter with a grid and
+    provides an interface to the solver. More complicated systems could be
+    modeled using a similar template. Peak-valued complex space vectors are
+    used. The space vector models are implemented in stationary coordinates.
+
+"""
+
+import numpy as np
+
+from gritulator._helpers import complex2abc
+from gritulator._utils import Bunch
+
+# %%
+class DCBusAndLFilterModel:
+    """
+    Continuous-time model for a stiff grid model with an RL impedance model.
+
+    Parameters
+    ----------
+    grid_filter : LFilter
+        RL line dynamic model.
+    grid_model : StiffSource
+        Constant voltage source model.
+    dc_model : DCBus
+        DC-bus voltage dynamics.
+    converter : Inverter | PWMInverter
+        Inverter model.
+
+    """
+    
+    def __init__(
+            self, grid_filter=None, grid_model=None,
+            dc_model=None, converter=None):
+        self.grid_filter = grid_filter
+        self.grid_model = grid_model
+        self.dc_model = dc_model
+        self.converter = converter
+
+        # Initial time
+        self.t0 = 0
+
+        # Store the solution in these lists
+        self.data = Bunch()
+        self.data.t, self.data.q = [], []
+        self.data.i_gs = []
+        self.data.u_dc = [] 
+        self.data.i_dc = []
+        
+    def get_initial_values(self):
+        """
+        Get the initial values.
+
+        Returns
+        -------
+        x0 : complex list, length 2
+            Initial values of the state variables.
+
+        """
+        x0 = [self.grid_filter.i_gs0, self.dc_model.u_dc0]
+        return x0
+
+    def set_initial_values(self, t0, x0):
+        """
+        Set the initial values.
+
+        Parameters
+        ----------
+        x0 : complex ndarray
+            Initial values of the state variables.
+
+        """
+        self.t0 = t0
+        self.grid_filter.i_gs0 = x0[0]
+        self.dc_model.u_dc0 = x0[1].real
+        self.converter.u_dc0 = x0[1].real
+        # calculation of converter-side voltage
+        u_cs0 = self.converter.ac_voltage(self.converter.q, x0[1].real)
+        # calculation of grid-side voltage
+        e_gs0 = self.grid_model.voltages(t0)
+        # update pcc voltage
+        self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(
+                                                x0[0], u_cs0, e_gs0)
+
+    def f(self, t, x):
+        """
+        Compute the complete state derivative list for the solver.
+
+        Parameters
+        ----------
+        t : float
+            Time.
+        x : complex ndarray
+            State vector.
+
+        Returns
+        -------
+        complex list
+            State derivatives.
+
+        """
+        # Unpack the states
+        i_gs, u_dc = x
+        # Interconnections: outputs for computing the state derivatives
+        u_cs = self.converter.ac_voltage(self.converter.q, u_dc)
+        e_gs = self.grid_model.voltages(t)
+        q = self.converter.q
+        i_g_abc = complex2abc(i_gs)
+        i_dc = self.dc_model.dc_current(i_g_abc, q) # DC-current
+        # State derivatives
+        rl_f = self.grid_filter.f(i_gs, u_cs, e_gs)
+        dc_f = self.dc_model.f(t, u_dc, i_dc)
+        # List of state derivatives
+        return [rl_f, dc_f]
+
+    def save(self, sol):
+        """
+        Save the solution.
+
+        Parameters
+        ----------
+        sol : Bunch object
+            Solution from the solver.
+
+        """
+        self.data.t.extend(sol.t)
+        self.data.i_gs.extend(sol.y[0])
+        self.data.u_dc.extend(sol.y[1].real)
+        self.data.q.extend(sol.q)
+        q_abc=complex2abc(np.asarray(sol.q))
+        i_c_abc=complex2abc(sol.y[0])
+        self.data.i_dc.extend(
+            q_abc[0]*i_c_abc[0] + q_abc[1]*i_c_abc[1] + q_abc[2]*i_c_abc[2])
+                                    
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format and post-process them.
+
+        """
+        # From lists to the ndarray
+        self.data.t = np.asarray(self.data.t)
+        self.data.i_gs = np.asarray(self.data.i_gs)
+        self.data.u_dc = np.asarray(self.data.u_dc)
+        self.data.q = np.asarray(self.data.q)
+        #self.data.theta = np.asarray(self.data.theta)
+        # Some useful variables
+        self.data.i_dc = np.asarray(self.data.i_dc)
+        self.data.e_gs = self.grid_model.voltages(self.data.t)
+        self.data.theta = np.mod(self.data.t*self.grid_model.w_N, 2*np.pi)
+        self.data.u_cs = self.converter.ac_voltage(
+            self.data.q, self.data.u_dc)
+        self.data.u_gs = self.grid_filter.pcc_voltages(
+            self.data.i_gs,
+            self.data.u_cs,
+            self.data.e_gs)
+
+
+# %%
+class DCBusAndLCLFilterModel:
+    """
+    Continuous-time model for a stiff grid model with an LCL impedance model.
+
+    Parameters
+    ----------
+    grid_filter : LCLFilter
+        LCL filter dynamic model.
+    grid_model : StiffSource
+        Constant voltage source model.
+    dc_model : DCBus
+        DC-bus voltage dynamics.
+    converter : Inverter | PWMInverter
+        Inverter model.
+
+    """
+    
+    def __init__(
+            self, grid_filter=None, grid_model=None, 
+            dc_model=None, converter=None):
+        self.grid_filter = grid_filter
+        self.grid_model = grid_model
+        self.dc_model = dc_model
+        self.converter = converter
+
+        # Initial time
+        self.t0 = 0
+
+        # Store the solution in these lists
+        self.data = Bunch()
+        self.data.t, self.data.q = [], []
+        self.data.i_gs = []
+        self.data.i_cs = []
+        self.data.u_fs = []
+        self.data.u_dc = [] 
+        self.data.i_dc = []
+        
+    def get_initial_values(self):
+        """
+        Get the initial values.
+
+        Returns
+        -------
+        x0 : complex list, length 4
+            Initial values of the state variables.
+
+        """
+        x0 = [
+            self.grid_filter.i_cs0,
+            self.grid_filter.u_fs0,
+            self.grid_filter.i_gs0,
+            self.dc_model.u_dc0]
+
+        return x0
+
+    def set_initial_values(self, t0, x0):
+        """
+        Set the initial values.
+
+        Parameters
+        ----------
+        x0 : complex ndarray
+            Initial values of the state variables.
+
+        """
+        self.t0 = t0
+        self.grid_filter.i_cs0 = x0[0]
+        self.grid_filter.u_fs0 = x0[1]
+        self.grid_filter.i_gs0 = x0[2]
+        self.dc_model.u_dc0 = x0[3].real
+        self.converter.u_dc0 = x0[3].real
+        # calculation of grid-side voltage
+        e_gs0 = self.grid_model.voltages(t0)
+        # update pcc voltage
+        self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(
+                                                    x0[2], x0[1], e_gs0)
+
+    def f(self, t, x):
+        """
+        Compute the complete state derivative list for the solver.
+
+        Parameters
+        ----------
+        t : float
+            Time.
+        x : complex ndarray
+            State vector.
+
+        Returns
+        -------
+        complex list
+            State derivatives.
+
+        """
+        # Unpack the states
+        i_cs, u_fs, i_gs, u_dc = x
+        # Interconnections: outputs for computing the state derivatives
+        u_cs = self.converter.ac_voltage(self.converter.q, self.converter.u_dc0)
+        e_gs = self.grid_model.voltages(t)
+        q = self.converter.q
+        i_c_abc = complex2abc(i_cs)
+        i_dc = self.dc_model.dc_current(i_c_abc, q) # DC-current
+        # State derivatives
+        lcl_f = self.grid_filter.f(i_cs, u_fs, i_gs, u_cs, e_gs)
+        dc_f = self.dc_model.f(t, u_dc, i_dc)
+        # List of state derivatives
+        all_f = [lcl_f[0], lcl_f[1], lcl_f[2], dc_f]
+        return all_f
+
+    def save(self, sol):
+        """
+        Save the solution.
+
+        Parameters
+        ----------
+        sol : Bunch object
+            Solution from the solver.
+
+        """
+        self.data.t.extend(sol.t)
+        self.data.i_cs.extend(sol.y[0])
+        self.data.u_fs.extend(sol.y[1])
+        self.data.i_gs.extend(sol.y[2])
+        self.data.u_dc.extend(sol.y[3].real)
+        self.data.q.extend(sol.q)
+        q_abc=complex2abc(np.asarray(sol.q))
+        i_c_abc=complex2abc(sol.y[0])
+        self.data.i_dc.extend(
+            q_abc[0]*i_c_abc[0] + q_abc[1]*i_c_abc[1] + q_abc[2]*i_c_abc[2])
+                                
+    def post_process(self):
+        """
+        Transform the lists to the ndarray format and post-process them.
+
+        """
+        # From lists to the ndarray
+        self.data.t = np.asarray(self.data.t)
+        self.data.i_cs = np.asarray(self.data.i_cs)
+        self.data.u_fs = np.asarray(self.data.u_fs)
+        self.data.i_gs = np.asarray(self.data.i_gs)
+        self.data.u_dc = np.asarray(self.data.u_dc)
+        self.data.q = np.asarray(self.data.q)
+        # Some useful variables
+        self.data.i_dc = np.asarray(self.data.i_dc)
+        self.data.e_gs = self.grid_model.voltages(self.data.t)
+        self.data.theta = np.mod(self.data.t*self.grid_model.w_N, 2*np.pi)
+        self.data.u_cs = self.converter.ac_voltage(self.data.q, self.data.u_dc)
+        self.data.u_gs = self.grid_filter.pcc_voltages(
+            self.data.i_gs,
+            self.data.u_fs,
+            self.data.e_gs)

From f2ce53042500c54d03921a03ef0591164059568e Mon Sep 17 00:00:00 2001
From: Mikko Saren <mikko.saren@aalto.fi>
Date: Wed, 12 Jun 2024 10:59:52 +0300
Subject: [PATCH 19/26] dc_bus propably ready

---
 motulator/grid/model/dc_bus/_dc_bus.py       | 23 ++++++++++++++------
 motulator/grid/model/dc_bus/_dc_dyn_model.py | 17 +++++++--------
 2 files changed, 24 insertions(+), 16 deletions(-)

diff --git a/motulator/grid/model/dc_bus/_dc_bus.py b/motulator/grid/model/dc_bus/_dc_bus.py
index 5417d2d17..195e38199 100644
--- a/motulator/grid/model/dc_bus/_dc_bus.py
+++ b/motulator/grid/model/dc_bus/_dc_bus.py
@@ -4,11 +4,14 @@
 A DC bus between an external current source or sink and a converter is modeled 
 considering an equivalent circuit comprising a capacitor and parallel resistor.
 """
+from types import SimpleNamespace
 
 from motulator.common.utils import complex2abc
 
+from motulator.common.model import Subsystem
+
 # %%
-class DCBus:
+class DCBus(Subsystem):
     """
     DC bus model
 
@@ -28,11 +31,17 @@ class DCBus:
 
     """
     def __init__(self, C_dc=1e-3, G_dc=0, i_ext=lambda t: 0, u_dc0 = 650):
-        self.C_dc = C_dc
-        self.G_dc = G_dc
-        self.i_ext = i_ext
+        super().__init__()
+        self.par = SimpleNamespace(C_dc=C_dc, G_dc=G_dc, i_ext=i_ext) 
+        self.udc0 = u_dc0
         # Initial values
-        self.u_dc0 = u_dc0
+        self.state = SimpleNamespace(u_dc = u_dc0)
+        self.sol_states = SimpleNamespace(u_dc = [])
+
+    def set_outputs(self, _):
+        """Set output variables."""
+        state, out = self.state, self.out
+        out.u_dc = state.u_dc
     
     @staticmethod
     def dc_current(i_c_abc, q):
@@ -80,7 +89,7 @@ def rhs(self, t, u_dc, i_dc):
 
         """
         # State derivative
-        du_dc = (self.i_ext(t) - i_dc - self.G_dc*u_dc)/self.C_dc
+        du_dc = (self.par.i_ext(t) - i_dc - self.par.G_dc*u_dc)/self.par.C_dc
         return du_dc
 
     def meas_dc_voltage(self):
@@ -93,4 +102,4 @@ def meas_dc_voltage(self):
             DC bus voltage (V)
     
         """
-        return self.u_dc0
+        return self.state.u_dc
diff --git a/motulator/grid/model/dc_bus/_dc_dyn_model.py b/motulator/grid/model/dc_bus/_dc_dyn_model.py
index bea8ca52d..1f65257b5 100644
--- a/motulator/grid/model/dc_bus/_dc_dyn_model.py
+++ b/motulator/grid/model/dc_bus/_dc_dyn_model.py
@@ -7,11 +7,10 @@
     used. The space vector models are implemented in stationary coordinates.
 
 """
-
+from types import SimpleNamespace
 import numpy as np
 
-from gritulator._helpers import complex2abc
-from gritulator._utils import Bunch
+from motulator.common.utils import complex2abc
 
 # %%
 class DCBusAndLFilterModel:
@@ -43,7 +42,7 @@ def __init__(
         self.t0 = 0
 
         # Store the solution in these lists
-        self.data = Bunch()
+        self.data = SimpleNamespace()
         self.data.t, self.data.q = [], []
         self.data.i_gs = []
         self.data.u_dc = [] 
@@ -84,7 +83,7 @@ def set_initial_values(self, t0, x0):
         self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(
                                                 x0[0], u_cs0, e_gs0)
 
-    def f(self, t, x):
+    def rhs(self, t, x):
         """
         Compute the complete state derivative list for the solver.
 
@@ -121,7 +120,7 @@ def save(self, sol):
 
         Parameters
         ----------
-        sol : Bunch object
+        sol : SimpleNamespace object
             Solution from the solver.
 
         """
@@ -187,7 +186,7 @@ def __init__(
         self.t0 = 0
 
         # Store the solution in these lists
-        self.data = Bunch()
+        self.data = SimpleNamespace()
         self.data.t, self.data.q = [], []
         self.data.i_gs = []
         self.data.i_cs = []
@@ -235,7 +234,7 @@ def set_initial_values(self, t0, x0):
         self.grid_filter.u_gs0 = self.grid_filter.pcc_voltages(
                                                     x0[2], x0[1], e_gs0)
 
-    def f(self, t, x):
+    def rhs(self, t, x):
         """
         Compute the complete state derivative list for the solver.
 
@@ -273,7 +272,7 @@ def save(self, sol):
 
         Parameters
         ----------
-        sol : Bunch object
+        sol : SimpleNamespace object
             Solution from the solver.
 
         """

From 4599ef3c0ed9bf4de4a8ced5e75206261433b652 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Thu, 13 Jun 2024 14:27:24 +0300
Subject: [PATCH 20/26] gfl example able to simulate, grid voltage still not
 saved

---
 .gitignore                                    |  2 ++
 ..._following_control_grid_converter_10kVA.py | 34 +++++++++++++++++--
 motulator/common/model/_simulation.py         |  1 +
 motulator/grid/control/_common.py             | 23 ++++---------
 .../control/grid_following/_grid_following.py | 13 +++++--
 motulator/grid/model/_grid_filter.py          | 21 ++++++++----
 motulator/grid/model/_grid_volt_source.py     |  8 ++---
 motulator/grid/utils/__init__.py              |  2 ++
 motulator/grid/utils/_plots.py                | 20 +++++------
 9 files changed, 82 insertions(+), 42 deletions(-)

diff --git a/.gitignore b/.gitignore
index 9ec9eca39..09f58ced8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,3 +14,5 @@ docs/source/auto_examples/
 docs/build/
 docs/source/autoapi/
 docs/source/sg_execution_times.rst
+
+*.json
\ No newline at end of file
diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index a25ca2461..50cfa9db5 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -14,11 +14,13 @@
 
 import time
 import numpy as np
+import matplotlib.pyplot as plt
 
 from motulator.grid import model
 import motulator.grid.control.grid_following as control
 #import motulator.grid.control.grid_following as control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
+from motulator.common.utils import complex2abc
 
 # To check the computation time of the program
 start_time = time.time()
@@ -41,7 +43,7 @@
 
 
 mdl = model.StiffSourceAndLFilterModel(converter, grid_filter, grid_model)
-mdl.pwm = None  # disable the PWM model
+
 
 
 
@@ -83,4 +85,32 @@
 # %%
 # Plot results in SI or per unit values.
 
-plot_grid(sim, base, plot_pcc_voltage=True)
+#plot_grid(sim, base, plot_pcc_voltage=True)
+
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))
+FS = 16 # Font size of the plots axis
+FL = 16 # Font size of the legends only
+LW = 3 # Line width in plots
+t_range = (0, 0.1)
+ctrl=sim.ctrl.data
+
+# Subplot 1: PCC voltage
+ax1.plot(ctrl.t, ctrl.u_g_abc/base.u, linewidth=LW)
+#ax1.plot(mdl.grid_filter.data.t, (mdl.grid_filter.data.e_gs)/base.u, linewidth=LW)
+ax1.legend([r'$u_g^a$',r'$u_g^b$',r'$u_g^c$'],
+            prop={'size': FL}, loc= 'upper right')
+ax1.set_xlim(t_range)
+ax1.set_xticklabels([])
+
+ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
+ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
+ax2.plot(ctrl.t, np.real(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+ax2.plot(ctrl.t, np.imag(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+#ax2.plot(mdl.t, mdl.iL, linewidth=LW) converter-side dc current for debug
+ax2.legend([r'$i_{c}^d$',r'$i_{c}^q$',r'$i_{c,ref}^d$',r'$i_{c,ref}^q$'],
+            prop={'size': FL}, loc= 'upper right')
+ax2.set_xlim(t_range)
+ax2.set_xticklabels([])
+
+plt.show()
+print(mdl.grid_filter.data.e_gs)
diff --git a/motulator/common/model/_simulation.py b/motulator/common/model/_simulation.py
index c2f575598..450b51f59 100644
--- a/motulator/common/model/_simulation.py
+++ b/motulator/common/model/_simulation.py
@@ -363,6 +363,7 @@ def save(self, sol):
                 for attr in vars(subsystem.sol_states):
                     subsystem.sol_states.__dict__[attr].extend(sol.y[index])
                     index += 1
+            
 
     def post_process_states(self):
         """Transform the lists to the ndarray format and post-process them."""
diff --git a/motulator/grid/control/_common.py b/motulator/grid/control/_common.py
index 3f3d32592..fc8d055f5 100644
--- a/motulator/grid/control/_common.py
+++ b/motulator/grid/control/_common.py
@@ -3,9 +3,9 @@
 from types import SimpleNamespace
 
 import numpy as np
-#from motulator.common.utils._utils import (abc2complex, complex2abc)
 #from motulator.common.control import (Ctrl, PICtrl)
 from motulator.common.utils import (complex2abc, abc2complex, wrap)
+from motulator.grid.utils import Bunch
 
 # %%
 class PWM:
@@ -597,7 +597,7 @@ class Ctrl:
     """Base class for the control system."""
 
     def __init__(self):
-        self.data = SimpleNamespace()  # Data store
+        self.data = Bunch()  # Data store
         self.clock = Clock()  # Digital clock
 
     def __call__(self, mdl):
@@ -616,7 +616,7 @@ def __call__(self, mdl):
         """
         raise NotImplementedError
 
-    def save(self, **kwargs):
+    def save(self, data):
         """
         Save the internal date of the control system.
 
@@ -626,14 +626,8 @@ def save(self, **kwargs):
             Contains the data to be saved.
 
         """
-        for name, arg in kwargs.items():
-            if not hasattr(self.data, name):
-                setattr(self.data, name, SimpleNamespace())
-            data = getattr(self.data, name)
-            for key, value in vars(arg).items():
-                if not hasattr(data, key):
-                    setattr(data, key, [])
-                getattr(data, key).append(value)
+        for key, value in data.items():
+            self.data.setdefault(key, []).extend([value])
 
     def post_process(self):
         """
@@ -643,11 +637,8 @@ def post_process(self):
         to simplify plotting and analysis of the stored data.
 
         """
-        self.data = SimpleNamespace()
-        for name, values in vars(self.data).items():
-            setattr(self.data, name, SimpleNamespace())
-            for key, value in vars(values).items():
-                setattr(getattr(self.data, name), key, np.asarray(value))
+        for key in self.data:
+            self.data[key] = np.asarray(self.data[key])
 
 
 # %%
diff --git a/motulator/grid/control/grid_following/_grid_following.py b/motulator/grid/control/grid_following/_grid_following.py
index 4f9e60419..5a063d7fa 100644
--- a/motulator/grid/control/grid_following/_grid_following.py
+++ b/motulator/grid/control/grid_following/_grid_following.py
@@ -5,8 +5,8 @@
 from collections.abc import Callable
 from dataclasses import dataclass, field
 
-from types import SimpleNamespace
-#from motulator._utils import Bunch
+#from types import SimpleNamespace
+from motulator.grid.utils import Bunch
 import numpy as np
 
 from motulator.common.utils._utils import abc2complex
@@ -202,7 +202,14 @@ def __call__(self, mdl):
                                            self.theta_p, self.w_g)
 
         # Data logging
-        self.save(fbk=mdl.fbk, ref=mdl.ref)
+        data = Bunch(
+            w_c = w_pll, theta_c = self.theta_p, u_c_ref = u_c_ref,
+            u_c = u_c, i_c = i_c, abs_u_g = abs_u_g,
+            d_abc_ref = d_abc_ref, i_c_ref = i_c_ref, u_dc = u_dc,
+            t = self.clock.t, p_g_ref = p_g_ref, u_dc_ref = u_dc_ref,
+            q_g_ref=q_g_ref, u_g = u_g, u_g_abc=u_g_abc,
+                     )
+        self.save(data)
 
         # Update the states
         self.theta_p = theta_pll
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index e23435fe8..a707d728f 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -39,12 +39,11 @@ class LFilter(Subsystem):
     def __init__(self, U_gN, L_f, R_f=0, L_g=0, R_g=0):
         super().__init__()
         self.par = SimpleNamespace(L_f=L_f, R_f=R_f, L_g=L_g, R_g=R_g)
-        #self.inp = SimpleNamespace(u_cs=0+0j)
+        self.inp = SimpleNamespace(u_cs=0+0j, e_gs=U_gN+0j)
         self.out = SimpleNamespace(u_gs=U_gN+0j)  # Needed for direct feedthrough
-        self.state = SimpleNamespace(i_gs=0)
+        self.state = SimpleNamespace(i_gs=0+0j)
         self.sol_states = SimpleNamespace(i_gs=[])
 
-
     def set_outputs(self, _):
         """Set output variables."""
         state, out = self.state, self.out
@@ -53,7 +52,6 @@ def set_outputs(self, _):
             self.state.i_gs)/(self.par.L_g+self.par.L_f)
         out.i_gs, out.i_cs, out.u_gs = state.i_gs, state.i_gs, u_gs
 
-
     def rhs(self):
         # pylint: disable=R0913
         """
@@ -69,8 +67,7 @@ def rhs(self):
         L_t = par.L_f + par.L_g
         R_t = par.R_f + par.R_g
         di_gs = (inp.u_cs - inp.e_gs - R_t*state.i_gs)/L_t
-        return di_gs
-
+        return [di_gs]
 
     def meas_currents(self):
         """
@@ -86,7 +83,6 @@ def meas_currents(self):
         i_g_abc = complex2abc(self.state.i_gs)
         return i_g_abc
 
-
     def meas_pcc_voltage(self):
         """
         Measure the phase voltages at PCC at the end of the sampling period.
@@ -101,6 +97,17 @@ def meas_pcc_voltage(self):
         u_g_abc = complex2abc(self.out.u_gs)
         return u_g_abc
 
+    def post_process_states(self):
+        """Post-process data."""
+        self.data.i_cs=self.data.i_gs
+
+    #def post_process_with_inputs(self):
+    #    """Post-process data with inputs."""
+    #    data=self.data
+    #    data.u_gs=(self.par.L_g*data.u_cs + self.par.L_f*data.e_gs +
+    #        (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
+    #        data.i_gs)/(self.par.L_g+self.par.L_f)
+
 
 # %%
 class LCLFilter:
diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index da47a6a75..f9eee7f20 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -35,13 +35,11 @@ def __init__(self, w_N=2*np.pi*50,
         super().__init__()
         self.par = SimpleNamespace(w_N=w_N, e_g_abs=e_g_abs)
 
-
     @property
     def w_N(self):
         """Grid constant frequency (rad/s)."""
         return self.par.w_N
 
-
     def voltages(self, t):
         """
         Compute the grid voltage in stationary frame.
@@ -68,12 +66,10 @@ def voltages(self, t):
         e_gs = abc2complex([e_g_a, e_g_b, e_g_c])
         return e_gs
 
-
     def set_outputs(self, t):
         """Set output variables."""
         self.out.e_gs = self.voltages(t)
 
-
     def meas_voltages(self, t):
         """
         Measure the phase voltages at the end of the sampling period.
@@ -93,6 +89,10 @@ def meas_voltages(self, t):
         e_g_abc = complex2abc(self.voltages(t))
         return e_g_abc
 
+    def post_process_states(self):
+        """Post-process the solution."""
+        self.data.e_gs=self.out.e_gs
+
 # %%
 class FlexSource(Subsystem):
     """
diff --git a/motulator/grid/utils/__init__.py b/motulator/grid/utils/__init__.py
index 5a76a4a5f..c16ea7042 100644
--- a/motulator/grid/utils/__init__.py
+++ b/motulator/grid/utils/__init__.py
@@ -2,6 +2,7 @@
 from motulator.grid.utils._helpers import BaseValues, NominalValues
 from motulator.grid.utils._plots import plot_grid
 from motulator.common.utils import Sequence, Step
+from motulator.grid.utils._utils import Bunch
 
 __all__ = [
     "BaseValues",
@@ -9,4 +10,5 @@
     "plot_grid",
     "Sequence",
     "Step",
+    "Bunch",
 ]
diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py
index 75259cb08..f13bbbaf9 100644
--- a/motulator/grid/utils/_plots.py
+++ b/motulator/grid/utils/_plots.py
@@ -46,7 +46,7 @@ def plot_grid(
     LW = 3 # Line width in plots
     
     
-    mdl = sim.mdl.data      # Continuous-time data
+    mdl = sim.mdl      # Continuous-time data
     ctrl = sim.ctrl.data    # Discrete-time data
     
     # Check if the time span was given
@@ -62,18 +62,18 @@ def plot_grid(
         pu_vals = True
         
     # 3-phase quantities
-    i_g_abc = complex2abc(mdl.i_gs).T
-    u_g_abc = complex2abc(mdl.u_gs).T
-    e_g_abc = complex2abc(mdl.e_gs).T
+    i_g_abc = complex2abc(mdl.grid_filter.data.i_gs).T
+    #u_g_abc = complex2abc(mdl.u_gs).T
+    #e_g_abc = complex2abc(mdl.e_gs).T
     
     # grid voltage magnitude
-    abs_e_g = np.abs(mdl.e_gs)
+    #abs_e_g = np.abs(mdl.e_gs)
     
     # Calculation of active and reactive powers
-    # p_g = 1.5*np.asarray(np.real(ctrl.u_g*np.conj(ctrl.i_c)))
-    # q_g = 1.5*np.asarray(np.imag(ctrl.u_g*np.conj(ctrl.i_c)))
-    p_g = 1.5*np.asarray(np.real(mdl.u_gs*np.conj(mdl.i_gs)))
-    q_g = 1.5*np.asarray(np.imag(mdl.u_gs*np.conj(mdl.i_gs)))
+    p_g = 1.5*np.asarray(np.real(ctrl.u_g*np.conj(ctrl.i_c)))
+    q_g = 1.5*np.asarray(np.imag(ctrl.u_g*np.conj(ctrl.i_c)))
+    #p_g = 1.5*np.asarray(np.real(mdl.u_gs*np.conj(mdl.i_gs)))
+    #q_g = 1.5*np.asarray(np.imag(mdl.u_gs*np.conj(mdl.i_gs)))
     p_g_ref = np.asarray(ctrl.p_g_ref)
     q_g_ref = np.asarray(ctrl.q_g_ref)
     
@@ -90,7 +90,7 @@ def plot_grid(
             ax1.set_xticklabels([])
         else:
             # Subplot 1: PCC voltage
-            ax1.plot(mdl.t, u_g_abc/base.u, linewidth=LW)
+            ax1.plot(ctrl.t, ctrl.u_g_abc/base.u, linewidth=LW)
             ax1.legend([r'$u_g^a$',r'$u_g^b$',r'$u_g^c$'],
                        prop={'size': FL}, loc= 'upper right')
             ax1.set_xlim(t_range)

From 4855bd42da6b0eb488ac9083c655f3821e6f1fb2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mikko=20Sar=C3=A9n?= <mikko.saren@aalto.fi>
Date: Mon, 17 Jun 2024 09:29:38 +0300
Subject: [PATCH 21/26] grid_volt_source, trying to set theta as state variable

---
 .DS_Store                                     | Bin 0 -> 6148 bytes
 ..._following_control_grid_converter_10kVA.py |  20 ++--
 motulator/.DS_Store                           | Bin 0 -> 6148 bytes
 motulator/grid/model/_grid_volt_source.py     |  26 ++++-
 motulator/grid/utils/_helpers.py              |   4 +
 motulator/grid/utils/_utils.py                |  93 ++++++++++++++++++
 6 files changed, 135 insertions(+), 8 deletions(-)
 create mode 100644 .DS_Store
 create mode 100644 motulator/.DS_Store
 create mode 100644 motulator/grid/utils/_utils.py

diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..2553326770162529a3e24a952428c43e1dc2f61c
GIT binary patch
literal 6148
zcmeHK-HOvd82zRjZQ`Q#LSb)&fVZV~cV%($CY!d%g09f4dZAL2rf4+Hl%%yvDfC+3
z!RPQvd>60u%+ErTt+!faKKSO$&vzy>XOfu=5ix%j?-MnM$c78Ft|8fAqFwkUtLYjC
z$YdXVIwhZ;(K*nVQNSp$?Fxu*_a=4dDV^b~UH^X1=n=JOOz~;Thf&HSL}Hs@cYv*5
z(&XeUiiy*Zib~wqN}Nq2`a~gmoKT8<HU4KMzD|#^i*%@mkp<K%_g;HjiFVL0rzy?Q
zcL?N^ml<>d#<#U`Z(v08VxnjIJ({FRF>SYhu-cAw<!aro+q?Gb;2oa_({MT~{4jqh
zO6Q!X(M3Lto+qREpmD3ivuT)Qqg)b_c!ZEwFOn?gbDz($xRBgHHQ0{r3>tfj#li7W
z(>;8!zihgT<AZk7y?=PLTsrpcyZ0XU2OrXL#=l?|FoykB+I5W=AXMKlvyVZZr92a=
z3FZx~!Do05o!(<DhES{}RU1=U){TiC%ewWkCMVRTC%{c_h*-&W28zbEPyzi8*}Gdf
zHZvrnfKlMzDj@a;7cS7#xKb#$4kYpl0JKmo4S8{w5FA&dr*WkaEik62Kt*Nxiop~e
z^{$TdG_DjXIx&6uVEV~S-%yx*I{J5II5AJ5NsR(VfwBTSs@M?cfA7!t|1!yxi~>f1
z|4IR7b%U;tDe1HI(&WTh>%!l_g^9S8LP<fQk7Jn;NAVt98rocu0eTu&3Q+?yKLSz)
LlNbg5ssg_O{FAGs

literal 0
HcmV?d00001

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 50cfa9db5..2163391b7 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -102,15 +102,21 @@
 ax1.set_xlim(t_range)
 ax1.set_xticklabels([])
 
-ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
-ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
-ax2.plot(ctrl.t, np.real(ctrl.i_c_ref/base.i), '--', linewidth=LW)
-ax2.plot(ctrl.t, np.imag(ctrl.i_c_ref/base.i), '--', linewidth=LW)
-#ax2.plot(mdl.t, mdl.iL, linewidth=LW) converter-side dc current for debug
-ax2.legend([r'$i_{c}^d$',r'$i_{c}^q$',r'$i_{c,ref}^d$',r'$i_{c,ref}^q$'],
-            prop={'size': FL}, loc= 'upper right')
+ # Subplot 2: Converter currents
+ax2.plot(ctrl.t, complex2abc(ctrl.u_gs).T/base.i, linewidth=LW)
+ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
+            ,prop={'size': FL}, loc= 'upper right')
 ax2.set_xlim(t_range)
 ax2.set_xticklabels([])
+# ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
+# ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
+# ax2.plot(ctrl.t, np.real(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+# ax2.plot(ctrl.t, np.imag(ctrl.i_c_ref/base.i), '--', linewidth=LW)
+# #ax2.plot(mdl.t, mdl.iL, linewidth=LW) converter-side dc current for debug
+# ax2.legend([r'$i_{c}^d$',r'$i_{c}^q$',r'$i_{c,ref}^d$',r'$i_{c,ref}^q$'],
+#             prop={'size': FL}, loc= 'upper right')
+# ax2.set_xlim(t_range)
+# ax2.set_xticklabels([])
 
 plt.show()
 print(mdl.grid_filter.data.e_gs)
diff --git a/motulator/.DS_Store b/motulator/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..26d2bd23f09eabbd8b3dd00475011b15548085d8
GIT binary patch
literal 6148
zcmeHKL2DC16n>MY>4qTWpwNpj;5F2=m16O-jrGv0D|%3g*>0$d+nJDT8YBevtUt!9
zzodVs{oc$}l14ohk@DX1=9_u*y;<g)WoCznRBxK?619maf-}~dD1Kwy&%S205NQCF
z8zZNZ`gnF|w8oae|EK`3UCXUmPKK7(@9RlnPt(GtVDRU01mC%Gh8j_WPRY_a9!usH
z_s8G3>1;9z<Ni!XMLsFJ-Ah&5sBbp5qOGVEeM~;vSyB$l>7+LpzvtE)Ym4-zK1ko@
z!&$%m@TE27AUDIYOUSb!C?DSCCbP4iotkXoaueSXHKS&~y*;1r9UgS_{`1{MN6!!U
zx*h#&|6s9bMvtC6eet^gr8qbCo0}mieAUh#1iXRosI1KYA{m>)nlr>XXOl7-L(K~s
zt}3G*qjy=wI-NOYc8kYvSmvQACJKlGqQIRg5O%kyb!Q{WK%#&sa9<SQ{XxSS<A9Y%
zw{@U!R{&rgVP~lG&m!mOfN{XeBYI%UrviPd3Reu}(-GH!7YD37`gBs^@}a_;Rk)!j
z_jZi09i3F{QA$xj6j)VY!#CS}{=fS3{eLw{TB3j`aIX|l_2cBYhbe`#b!~Ed)@t}C
qI2-dSkH1n-aYr$@d=ziPoguEd0*nJz9ua}bkARaw3Q^#nDsT<PmwW~Q

literal 0
HcmV?d00001

diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index f9eee7f20..a474b3837 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -7,6 +7,11 @@
 synchronous generator are considered. In this module, all space vectors are in
 stationary coordinates.
 
+The grid angle theta_g is used 
+    as state variables.
+
+        theta_g = w_N * t
+
 """
 from types import SimpleNamespace
 
@@ -34,12 +39,31 @@ def __init__(self, w_N=2*np.pi*50,
                  e_g_abs=lambda t: 400*np.sqrt(2/3)):
         super().__init__()
         self.par = SimpleNamespace(w_N=w_N, e_g_abs=e_g_abs)
+        # states
+        self.state = SimpleNamespace(theta_g=0)
+        # Store the solutions in these lists
+        self.sol_states = SimpleNamespace(theta_g=[]) 
 
     @property
     def w_N(self):
         """Grid constant frequency (rad/s)."""
+        if callable(self.par.w_N):
+            return self.par.w_N()
         return self.par.w_N
-
+    
+    def rhs(self):
+        """
+        Compute the state derivatives.
+        
+        Returns
+        -------
+        list, length 1
+            Time derivatives of the state vector.
+            
+        """
+        dtheta_g = self.par.w_N
+        return [dtheta_g]
+        
     def voltages(self, t):
         """
         Compute the grid voltage in stationary frame.
diff --git a/motulator/grid/utils/_helpers.py b/motulator/grid/utils/_helpers.py
index c7443e6a5..7da6e7fa2 100644
--- a/motulator/grid/utils/_helpers.py
+++ b/motulator/grid/utils/_helpers.py
@@ -1,7 +1,11 @@
 """Helper functions and classes."""
 
 # %%
+from abc import ABC
 from dataclasses import dataclass
+# Note: Union can be replaced by "|" in Python 3.10
+from typing import Callable, Union
+
 import numpy as np
 
 # %%
diff --git a/motulator/grid/utils/_utils.py b/motulator/grid/utils/_utils.py
new file mode 100644
index 000000000..093621172
--- /dev/null
+++ b/motulator/grid/utils/_utils.py
@@ -0,0 +1,93 @@
+"""Utils."""
+
+
+# %%
+# The Bunch class is from sklearn.utils:
+#
+#    https://github.com/scikit-learn/scikit-learn
+#
+# Its license is provided below.
+#
+# ----------------------------------------------------------------------------
+# BSD 3-Clause License
+#
+# Copyright (c) 2007-2021 The scikit-learn developers.
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of the copyright holder nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+# ----------------------------------------------------------------------------
+class Bunch(dict):
+    """
+    Container object exposing keys as attributes.
+
+    Bunch objects are sometimes used as an output for functions and methods.
+    They extend dictionaries by enabling values to be accessed by key,
+    `bunch["value_key"]`, or by an attribute, `bunch.value_key`.
+
+    Examples
+    --------
+    >>> from sklearn.utils import Bunch
+    >>> b = Bunch(a=1, b=2)
+    >>> b['b']
+    2
+    >>> b.b
+    2
+    >>> b.a = 3
+    >>> b['a']
+    3
+    >>> b.c = 6
+    >>> b['c']
+    6
+
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(kwargs)
+
+    def __setattr__(self, key, value):
+        self[key] = value
+
+    def __dir__(self):
+        return self.keys()
+
+    def __getattr__(self, key):
+        try:
+            return self[key]
+        except KeyError:
+            # pylint: disable=raise-missing-from
+            raise AttributeError(key)
+
+    def __setstate__(self, state):
+        # Bunch pickles generated with scikit-learn 0.16.* have an non
+        # empty __dict__. This causes a surprising behaviour when
+        # loading these pickles scikit-learn 0.17: reading bunch.key
+        # uses __dict__ but assigning to bunch.key use __setattr__ and
+        # only changes bunch['key']. More details can be found at:
+        # https://github.com/scikit-learn/scikit-learn/issues/6196.
+        # Overriding __setstate__ to be a noop has the effect of
+        # ignoring the pickled __dict__
+        pass

From c5bd1d9b3359f4dd2311de5502ee306463bc6875 Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 17 Jun 2024 10:30:30 +0300
Subject: [PATCH 22/26] grid voltage angle as state variable, pcc voltage
 calculation

---
 ..._following_control_grid_converter_10kVA.py |  6 ++---
 motulator/common/model/_simulation.py         |  2 +-
 motulator/grid/model/_grid_filter.py          | 12 ++++-----
 motulator/grid/model/_grid_volt_source.py     | 27 ++++++++++---------
 4 files changed, 24 insertions(+), 23 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 2163391b7..30a428d06 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -103,9 +103,9 @@
 ax1.set_xticklabels([])
 
  # Subplot 2: Converter currents
-ax2.plot(ctrl.t, complex2abc(ctrl.u_gs).T/base.i, linewidth=LW)
-ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
-            ,prop={'size': FL}, loc= 'upper right')
+ax2.plot(mdl.grid_filter.data.t, mdl.grid_filter.data.e_gs/base.u, linewidth=LW)
+#ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
+#            ,prop={'size': FL}, loc= 'upper right')
 ax2.set_xlim(t_range)
 ax2.set_xticklabels([])
 # ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
diff --git a/motulator/common/model/_simulation.py b/motulator/common/model/_simulation.py
index 450b51f59..a798418b2 100644
--- a/motulator/common/model/_simulation.py
+++ b/motulator/common/model/_simulation.py
@@ -363,7 +363,7 @@ def save(self, sol):
                 for attr in vars(subsystem.sol_states):
                     subsystem.sol_states.__dict__[attr].extend(sol.y[index])
                     index += 1
-            
+
 
     def post_process_states(self):
         """Transform the lists to the ndarray format and post-process them."""
diff --git a/motulator/grid/model/_grid_filter.py b/motulator/grid/model/_grid_filter.py
index a707d728f..74a1845a4 100644
--- a/motulator/grid/model/_grid_filter.py
+++ b/motulator/grid/model/_grid_filter.py
@@ -101,12 +101,12 @@ def post_process_states(self):
         """Post-process data."""
         self.data.i_cs=self.data.i_gs
 
-    #def post_process_with_inputs(self):
-    #    """Post-process data with inputs."""
-    #    data=self.data
-    #    data.u_gs=(self.par.L_g*data.u_cs + self.par.L_f*data.e_gs +
-    #        (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
-    #        data.i_gs)/(self.par.L_g+self.par.L_f)
+    def post_process_with_inputs(self):
+        """Post-process data with inputs."""
+        data=self.data
+        data.u_gs=(self.par.L_g*data.u_cs + self.par.L_f*data.e_gs +
+            (self.par.R_g*self.par.L_f - self.par.R_f*self.par.L_g)*
+            data.i_gs)/(self.par.L_g+self.par.L_f)
 
 
 # %%
diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index a474b3837..086b4846b 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -42,7 +42,7 @@ def __init__(self, w_N=2*np.pi*50,
         # states
         self.state = SimpleNamespace(theta_g=0)
         # Store the solutions in these lists
-        self.sol_states = SimpleNamespace(theta_g=[]) 
+        self.sol_states = SimpleNamespace(theta_g=[])
 
     @property
     def w_N(self):
@@ -50,7 +50,7 @@ def w_N(self):
         if callable(self.par.w_N):
             return self.par.w_N()
         return self.par.w_N
-    
+
     def rhs(self):
         """
         Compute the state derivatives.
@@ -63,8 +63,8 @@ def rhs(self):
         """
         dtheta_g = self.par.w_N
         return [dtheta_g]
-        
-    def voltages(self, t):
+
+    def voltages(self, t, theta_g):
         """
         Compute the grid voltage in stationary frame.
            
@@ -72,6 +72,8 @@ def voltages(self, t):
         ----------
         t : float
             Time (s).
+        theta_g : float
+            Grid voltage angle (rad)
 
         Returns
         -------
@@ -79,20 +81,18 @@ def voltages(self, t):
             grid complex voltage (V).
 
         """
-        # Calculation of theta angle based on time and grid fixed frequency
-        theta = self.w_N*t
 
-        # Calculation of the three-phase voltage
-        e_g_a = self.par.e_g_abs(t)*np.cos(theta)
-        e_g_b = self.par.e_g_abs(t)*np.cos(theta-2*np.pi/3)
-        e_g_c = self.par.e_g_abs(t)*np.cos(theta-4*np.pi/3)
+        # Calculation of the three-phase voltages
+        e_g_a = self.par.e_g_abs(t)*np.cos(theta_g)
+        e_g_b = self.par.e_g_abs(t)*np.cos(theta_g-2*np.pi/3)
+        e_g_c = self.par.e_g_abs(t)*np.cos(theta_g-4*np.pi/3)
 
         e_gs = abc2complex([e_g_a, e_g_b, e_g_c])
         return e_gs
 
     def set_outputs(self, t):
         """Set output variables."""
-        self.out.e_gs = self.voltages(t)
+        self.out.e_gs = self.voltages(t, self.state.theta_g)
 
     def meas_voltages(self, t):
         """
@@ -110,12 +110,13 @@ def meas_voltages(self, t):
 
         """
         # Grid voltage
-        e_g_abc = complex2abc(self.voltages(t))
+        e_g_abc = complex2abc(self.voltages(t, self.state.theta_g))
         return e_g_abc
 
     def post_process_states(self):
         """Post-process the solution."""
-        self.data.e_gs=self.out.e_gs
+        self.data.e_gs=self.voltages(self.data.t, self.data.theta_g)
+
 
 # %%
 class FlexSource(Subsystem):

From d62775654550b3c5ba5debf145401415e3ae66db Mon Sep 17 00:00:00 2001
From: maattaj11 <juho.k.maatta@aalto.fi>
Date: Mon, 17 Jun 2024 14:17:29 +0300
Subject: [PATCH 23/26] changed grid angle to use exponential form, working
 grid plots

---
 ..._following_control_grid_converter_10kVA.py | 45 +------------
 motulator/grid/model/_grid_volt_source.py     | 13 ++--
 motulator/grid/utils/_plots.py                | 67 +++++++++----------
 3 files changed, 39 insertions(+), 86 deletions(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 30a428d06..710c95e3f 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -14,13 +14,11 @@
 
 import time
 import numpy as np
-import matplotlib.pyplot as plt
 
 from motulator.grid import model
 import motulator.grid.control.grid_following as control
 #import motulator.grid.control.grid_following as control
 from motulator.grid.utils import BaseValues, NominalValues, plot_grid
-from motulator.common.utils import complex2abc
 
 # To check the computation time of the program
 start_time = time.time()
@@ -41,12 +39,9 @@
 grid_model = model.StiffSource(w_N=2*np.pi*50)
 converter = model.Inverter(u_dc=650)
 
-
 mdl = model.StiffSourceAndLFilterModel(converter, grid_filter, grid_model)
 
 
-
-
 # %%
 # Configure the control system.
 
@@ -81,42 +76,4 @@
 # Print the execution time
 print('\nExecution time: {:.2f} s'.format((time.time() - start_time)))
 
-
-# %%
-# Plot results in SI or per unit values.
-
-#plot_grid(sim, base, plot_pcc_voltage=True)
-
-fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))
-FS = 16 # Font size of the plots axis
-FL = 16 # Font size of the legends only
-LW = 3 # Line width in plots
-t_range = (0, 0.1)
-ctrl=sim.ctrl.data
-
-# Subplot 1: PCC voltage
-ax1.plot(ctrl.t, ctrl.u_g_abc/base.u, linewidth=LW)
-#ax1.plot(mdl.grid_filter.data.t, (mdl.grid_filter.data.e_gs)/base.u, linewidth=LW)
-ax1.legend([r'$u_g^a$',r'$u_g^b$',r'$u_g^c$'],
-            prop={'size': FL}, loc= 'upper right')
-ax1.set_xlim(t_range)
-ax1.set_xticklabels([])
-
- # Subplot 2: Converter currents
-ax2.plot(mdl.grid_filter.data.t, mdl.grid_filter.data.e_gs/base.u, linewidth=LW)
-#ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
-#            ,prop={'size': FL}, loc= 'upper right')
-ax2.set_xlim(t_range)
-ax2.set_xticklabels([])
-# ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
-# ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
-# ax2.plot(ctrl.t, np.real(ctrl.i_c_ref/base.i), '--', linewidth=LW)
-# ax2.plot(ctrl.t, np.imag(ctrl.i_c_ref/base.i), '--', linewidth=LW)
-# #ax2.plot(mdl.t, mdl.iL, linewidth=LW) converter-side dc current for debug
-# ax2.legend([r'$i_{c}^d$',r'$i_{c}^q$',r'$i_{c,ref}^d$',r'$i_{c,ref}^q$'],
-#             prop={'size': FL}, loc= 'upper right')
-# ax2.set_xlim(t_range)
-# ax2.set_xticklabels([])
-
-plt.show()
-print(mdl.grid_filter.data.e_gs)
+plot_grid(sim=sim, base=base, plot_pcc_voltage=False)
diff --git a/motulator/grid/model/_grid_volt_source.py b/motulator/grid/model/_grid_volt_source.py
index 086b4846b..07310e828 100644
--- a/motulator/grid/model/_grid_volt_source.py
+++ b/motulator/grid/model/_grid_volt_source.py
@@ -40,9 +40,9 @@ def __init__(self, w_N=2*np.pi*50,
         super().__init__()
         self.par = SimpleNamespace(w_N=w_N, e_g_abs=e_g_abs)
         # states
-        self.state = SimpleNamespace(theta_g=0)
+        self.state = SimpleNamespace(exp_j_theta_g=complex(1))
         # Store the solutions in these lists
-        self.sol_states = SimpleNamespace(theta_g=[])
+        self.sol_states = SimpleNamespace(exp_j_theta_g=[])
 
     @property
     def w_N(self):
@@ -61,8 +61,8 @@ def rhs(self):
             Time derivatives of the state vector.
             
         """
-        dtheta_g = self.par.w_N
-        return [dtheta_g]
+        d_exp_j_theta_g = 1j*self.par.w_N*self.state.exp_j_theta_g
+        return [d_exp_j_theta_g]
 
     def voltages(self, t, theta_g):
         """
@@ -92,7 +92,7 @@ def voltages(self, t, theta_g):
 
     def set_outputs(self, t):
         """Set output variables."""
-        self.out.e_gs = self.voltages(t, self.state.theta_g)
+        self.out.e_gs = self.voltages(t, np.angle(self.state.exp_j_theta_g))
 
     def meas_voltages(self, t):
         """
@@ -110,11 +110,12 @@ def meas_voltages(self, t):
 
         """
         # Grid voltage
-        e_g_abc = complex2abc(self.voltages(t, self.state.theta_g))
+        e_g_abc = complex2abc(self.voltages(t, np.angle(self.state.exp_j_theta_g)))
         return e_g_abc
 
     def post_process_states(self):
         """Post-process the solution."""
+        self.data.theta_g = np.angle(self.data.exp_j_theta_g)
         self.data.e_gs=self.voltages(self.data.t, self.data.theta_g)
 
 
diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py
index f13bbbaf9..3ebcdc395 100644
--- a/motulator/grid/utils/_plots.py
+++ b/motulator/grid/utils/_plots.py
@@ -5,7 +5,7 @@
 import matplotlib.pyplot as plt
 from cycler import cycler
 
-from motulator.common.utils import complex2abc
+from motulator.common.utils import (complex2abc, wrap)
 from types import SimpleNamespace
 
 # Plotting parameters
@@ -44,14 +44,14 @@ def plot_grid(
     FS = 16 # Font size of the plots axis
     FL = 16 # Font size of the legends only
     LW = 3 # Line width in plots
-    
-    
+
+
     mdl = sim.mdl      # Continuous-time data
     ctrl = sim.ctrl.data    # Discrete-time data
-    
+
     # Check if the time span was given
     if t_range is None:
-        t_range = (0, mdl.t[-1]) # Time span
+        t_range = (0, mdl.converter.data.t[-1]) # Time span
 
     # Check if the base values were given
     if base is None:
@@ -60,30 +60,26 @@ def plot_grid(
         base = SimpleNamespace(w=1, u=1, i=1, p=1000)
     else:
         pu_vals = True
-        
+
     # 3-phase quantities
-    i_g_abc = complex2abc(mdl.grid_filter.data.i_gs).T
-    #u_g_abc = complex2abc(mdl.u_gs).T
-    #e_g_abc = complex2abc(mdl.e_gs).T
-    
-    # grid voltage magnitude
-    #abs_e_g = np.abs(mdl.e_gs)
-    
+    i_c_abc = complex2abc(mdl.converter.data.i_cs).T
+    e_g_abc = complex2abc(mdl.grid_model.data.e_gs).T
+
     # Calculation of active and reactive powers
-    p_g = 1.5*np.asarray(np.real(ctrl.u_g*np.conj(ctrl.i_c)))
-    q_g = 1.5*np.asarray(np.imag(ctrl.u_g*np.conj(ctrl.i_c)))
-    #p_g = 1.5*np.asarray(np.real(mdl.u_gs*np.conj(mdl.i_gs)))
-    #q_g = 1.5*np.asarray(np.imag(mdl.u_gs*np.conj(mdl.i_gs)))
+    #p_g = 1.5*np.asarray(np.real(ctrl.u_g*np.conj(ctrl.i_c)))
+    #q_g = 1.5*np.asarray(np.imag(ctrl.u_g*np.conj(ctrl.i_c)))
+    p_g = 1.5*np.asarray(np.real(mdl.grid_filter.data.u_gs*np.conj(mdl.grid_filter.data.i_gs)))
+    q_g = 1.5*np.asarray(np.imag(mdl.grid_filter.data.u_gs*np.conj(mdl.grid_filter.data.i_gs)))
     p_g_ref = np.asarray(ctrl.p_g_ref)
     q_g_ref = np.asarray(ctrl.q_g_ref)
-    
+
     # %%
     fig, (ax1, ax2,ax3) = plt.subplots(3, 1, figsize=(10, 7))
 
     if sim.ctrl.on_v_dc==False:
         if plot_pcc_voltage == False:
             # Subplot 1: Grid voltage
-            ax1.plot(mdl.t, e_g_abc/base.u, linewidth=LW)
+            ax1.plot(mdl.grid_model.data.t, e_g_abc/base.u, linewidth=LW)
             ax1.legend([r'$e_g^a$',r'$e_g^b$',r'$e_g^c$'],
                        prop={'size': FL}, loc= 'upper right')
             ax1.set_xlim(t_range)
@@ -98,31 +94,31 @@ def plot_grid(
             #ax1.set_ylabel('Grid voltage (V)')
     else:
         # Subplot 1: DC-bus voltage
-        ax1.plot(mdl.t, mdl.u_dc/base.u, linewidth=LW)
+        ax1.plot(mdl.converter.data.t, mdl.converter.data.u_dc/base.u, linewidth=LW)
         ax1.plot(ctrl.t, ctrl.u_dc_ref/base.u, '--', linewidth=LW)
         ax1.legend([r'$u_{dc}$',r'$u_{dc,ref}$'],
                    prop={'size': FL}, loc= 'upper right')
         ax1.set_xlim(t_range)
         ax1.set_xticklabels([])
         #ax1.set_ylabel('DC-bus voltage (V)')
-    
+
     # Subplot 2: Converter currents
-    ax2.plot(mdl.t, i_g_abc/base.i, linewidth=LW)
-    ax2.legend([r'$i_g^a$',r'$i_g^b$',r'$i_g^c$']
+    ax2.plot(mdl.converter.data.t, i_c_abc/base.i, linewidth=LW)
+    ax2.legend([r'$i_c^a$',r'$i_c^b$',r'$i_c^c$']
                ,prop={'size': FL}, loc= 'upper right')
     ax2.set_xlim(t_range)
     ax2.set_xticklabels([])
-    
+
     if plot_w:
         # Subplot 3: Grid and converter frequencies
-        ax3.plot(mdl.t, mdl.w_g/base.w, linewidth=LW)
+        ax3.plot(mdl.grid_model.data.t, mdl.grid_model.par.w_N/base.w, linewidth=LW)
         ax3.plot(ctrl.t, ctrl.w_c/base.w, '--', linewidth=LW)
         ax3.legend([r'$\omega_{g}$',r'$\omega_{c}$']
                    ,prop={'size': FL}, loc= 'upper right')
         ax3.set_xlim(t_range)
     else:
         # Subplot 3: Phase angles
-        ax3.plot(mdl.t, mdl.theta, linewidth=LW)
+        ax3.plot(mdl.grid_model.data.t, mdl.grid_model.data.theta_g, linewidth=LW)
         ax3.plot(ctrl.t, ctrl.theta_c, '--', linewidth=LW)
         ax3.legend([r'$\theta_{g}$',r'$\theta_{c}$']
                    ,prop={'size': FL}, loc= 'upper right')
@@ -162,21 +158,21 @@ def plot_grid(
     ax3.grid()
     #plt.show()
 
-    
+
     # %%
     # Second figure
     fig, (ax1, ax2,ax3) = plt.subplots(3, 1, figsize=(10, 7))
 
     # Subplot 1: Active and reactive power
-    ax1.plot(mdl.t, p_g/base.p, linewidth=LW)
-    ax1.plot(mdl.t, q_g/base.p, linewidth=LW)
+    ax1.plot(mdl.grid_filter.data.t, p_g/base.p, linewidth=LW)
+    ax1.plot(mdl.grid_filter.data.t, q_g/base.p, linewidth=LW)
     ax1.plot(ctrl.t, (p_g_ref/base.p), '--', linewidth=LW)
     ax1.plot(ctrl.t, (q_g_ref/base.p), '--', linewidth=LW)
     ax1.legend([r'$p_{g}$',r'$q_{g}$',r'$p_{g,ref}$',r'$q_{g,ref}$'],
                prop={'size': FL}, loc= 'upper right')
     ax1.set_xlim(t_range)
     ax1.set_xticklabels([])
-    
+
     # Subplot 2: Converter currents
     ax2.plot(ctrl.t, np.real(ctrl.i_c/base.i), linewidth=LW)
     ax2.plot(ctrl.t, np.imag(ctrl.i_c/base.i), linewidth=LW)
@@ -189,14 +185,14 @@ def plot_grid(
     ax2.set_xticklabels([])
 
     # Subplot 3: Converter voltage reference and grid voltage
-    ax3.plot(ctrl.t,np.real(ctrl.u_c/base.u), 
+    ax3.plot(ctrl.t,np.real(ctrl.u_c/base.u),
             ctrl.t,np.imag(ctrl.u_c/base.u), linewidth=LW)
-    ax3.plot(mdl.t,np.real(abs_e_g/base.u),'k--', 
+    ax3.plot(mdl.grid_model.data.t, np.absolute(mdl.grid_model.data.e_gs/base.u),'k--',
              linewidth=LW)
-    ax3.legend([r'$u_{c}^d$', r'$u_{c}^q$', r'$|e_{g}|$'], 
+    ax3.legend([r'$u_{c}^d$', r'$u_{c}^q$', r'$|e_{g}|$'],
                 prop={'size': FS}, loc= 'upper right')
     ax3.set_xlim(t_range)
-    
+
     # Change font size
     for item in ([ax2.title, ax2.xaxis.label, ax2.yaxis.label] +
              ax2.get_xticklabels() + ax2.get_yticklabels()):
@@ -225,8 +221,7 @@ def plot_grid(
     plt.tight_layout()
     plt.grid()
     ax3.grid()
-    plt.show()    
-
+    plt.show()
 
 
 # %%

From 8111f0d43a380dd0f6d3aa9796cfd7db03078bd6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mikko=20Sar=C3=A9n?= <mikko.saren@aalto.fi>
Date: Tue, 18 Jun 2024 09:44:16 +0300
Subject: [PATCH 24/26] added post process to dc_bus

---
 motulator/common/control/_control.py   | 2 +-
 motulator/grid/model/dc_bus/_dc_bus.py | 6 ++++++
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/motulator/common/control/_control.py b/motulator/common/control/_control.py
index 23e031cfd..95bedf777 100644
--- a/motulator/common/control/_control.py
+++ b/motulator/common/control/_control.py
@@ -381,7 +381,7 @@ def update(self, T_s, u, w):
 
 
 # %%
-class ComplexFFPICtrl:
+class ComplexFFPIController:
     """
     2DOF Synchronous-frame complex-vector PI controller with feedforward.
 
diff --git a/motulator/grid/model/dc_bus/_dc_bus.py b/motulator/grid/model/dc_bus/_dc_bus.py
index 195e38199..7106ec6b1 100644
--- a/motulator/grid/model/dc_bus/_dc_bus.py
+++ b/motulator/grid/model/dc_bus/_dc_bus.py
@@ -103,3 +103,9 @@ def meas_dc_voltage(self):
     
         """
         return self.state.u_dc
+
+    def post_process_states(self):
+        """Post-process the solution."""
+        data, par = self.data, self.par
+        data.u_dc = self.state.u_dc
+        data.i_dc = par.i_ext(self.data.t) - par.G_dc*self.data.u_dc
\ No newline at end of file

From b1ef1b67fcf82ad40ca3475b4dd695565ab31413 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mikko=20Sar=C3=A9n?= <mikko.saren@aalto.fi>
Date: Tue, 18 Jun 2024 13:05:39 +0300
Subject: [PATCH 25/26] adding set inputs function

---
 motulator/grid/model/dc_bus/_dc_bus.py | 45 +++++++-------------------
 1 file changed, 11 insertions(+), 34 deletions(-)

diff --git a/motulator/grid/model/dc_bus/_dc_bus.py b/motulator/grid/model/dc_bus/_dc_bus.py
index 7106ec6b1..9d2855c14 100644
--- a/motulator/grid/model/dc_bus/_dc_bus.py
+++ b/motulator/grid/model/dc_bus/_dc_bus.py
@@ -32,45 +32,24 @@ class DCBus(Subsystem):
     """
     def __init__(self, C_dc=1e-3, G_dc=0, i_ext=lambda t: 0, u_dc0 = 650):
         super().__init__()
-        self.par = SimpleNamespace(C_dc=C_dc, G_dc=G_dc, i_ext=i_ext) 
+        self.par = SimpleNamespace(C_dc=C_dc, G_dc=G_dc, i_ext=i_ext)
         self.udc0 = u_dc0
+        self.i_ext = i_ext
         # Initial values
         self.state = SimpleNamespace(u_dc = u_dc0)
+        self.inp = SimpleNamespace(i_dc = 0, i_ext = i_ext(0))
         self.sol_states = SimpleNamespace(u_dc = [])
 
-    def set_outputs(self, _):
+    def set_outputs(self, t):
         """Set output variables."""
         state, out = self.state, self.out
         out.u_dc = state.u_dc
+
+    def set_inputs(self, t):
+        """Set input variables."""
+        self.inp.i_ext = self.i_ext(t)
     
-    @staticmethod
-    def dc_current(i_c_abc, q):
-        """
-        Compute the converter DC current from the switching states and phase 
-        currents.
-    
-        Parameters
-        ----------
-        i_c_abc : ndarray, shape (3,)
-            Phase currents (A).
-        q : complex ndarray, shape (3,)
-            Switching state vectors corresponding to the switching instants.
-            For example, the switching state q[1] is applied at the interval
-            t_n_sw[1].
-    
-        Returns
-        -------
-        i_dc: float
-            Converter DC current (A)
-    
-        """
-        # Duty ratio back into three-phase ratios
-        q_abc = complex2abc(q)
-        # Dot product
-        i_dc = q_abc[0]*i_c_abc[0] + q_abc[1]*i_c_abc[1] + q_abc[2]*i_c_abc[2]
-        return i_dc
-    
-    def rhs(self, t, u_dc, i_dc):
+    def rhs(self):
         """
         Compute the state derivatives.
 
@@ -89,7 +68,7 @@ def rhs(self, t, u_dc, i_dc):
 
         """
         # State derivative
-        du_dc = (self.par.i_ext(t) - i_dc - self.par.G_dc*u_dc)/self.par.C_dc
+        du_dc = (self.inp.i_ext - self.inp.i_dc - self.par.G_dc*self.state.u_dc)/self.par.C_dc
         return du_dc
 
     def meas_dc_voltage(self):
@@ -106,6 +85,4 @@ def meas_dc_voltage(self):
 
     def post_process_states(self):
         """Post-process the solution."""
-        data, par = self.data, self.par
-        data.u_dc = self.state.u_dc
-        data.i_dc = par.i_ext(self.data.t) - par.G_dc*self.data.u_dc
\ No newline at end of file
+        self.data.u_dc = self.state.u_dc

From 9a1b1bbb8556310e2079511b59054bdaf3e2ccbd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mikko=20Sar=C3=A9n?= <mikko.saren@aalto.fi>
Date: Tue, 18 Jun 2024 13:38:57 +0300
Subject: [PATCH 26/26] added placeholder for the PWM parameter in example

---
 .../grid/plot_grid_following_control_grid_converter_10kVA.py    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
index 710c95e3f..40901a29b 100644
--- a/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
+++ b/examples/grid/plot_grid_following_control_grid_converter_10kVA.py
@@ -66,10 +66,10 @@
 e_g_abs_var =  lambda t: np.sqrt(2/3)*400
 mdl.grid_model.e_g_abs = e_g_abs_var # grid voltage magnitude
 
-
 # %%
 # Create the simulation object and simulate it.
 
+#mdl.pwm = model.CarrierComparison()  # Enable the PWM model
 sim = model.Simulation(mdl, ctrl)
 sim.simulate(t_stop = .1)