PyIndMach012: an example of user-model using DSS-Python

This example runs a modified example from the OpenDSS distribution for the induction machine model with a sample PyIndMach012 implementation, written in Python, and the original, built-in IndMach012.

Check the file for more comments. Comparing it to the Pascal code for IndMach012 can be useful to understand some of the inner workings of OpenDSS.

The user-model code in DSS-Python was supposed to grow other features, but the effort hasn’t been continued for many reasons. The code for generator user-models has been stable for many releases. We believe it can be used to develop new ideas before committing the final model in a traditional DLL user-model.

The outputs for this notebook should be listed below, but you can also open and then run this notebook on Google Colab for a quick overview if you don’t want to set up a local environment: Open in Colab. The next cell installs DSS-Python and dependencies, and downloads the extra files.

[ ]:
import os, subprocess
if os.getenv("COLAB_RELEASE_TAG"):
    import urllib.request
    print(subprocess.check_output('pip install dss-python[plot]', shell=True).decode())
    urllib.request.urlretrieve("", "")
    urllib.request.urlretrieve("", "master.dss")

import os
import numpy as np
from matplotlib import pyplot as plt
from dss.UserModels import GenUserModel # used to get the DLL path
import PyIndMach012 # we need to import the model so it gets registered

The model class

Type:        module
String form: <module 'PyIndMach012' from '/home/meira/Projects/dss/dss_python/docs/examples/UserModels/PyIndMach012/'>
File:        ~/Projects/dss/dss_python/docs/examples/UserModels/PyIndMach012/
A `dss_python` User-Model implementation of the IndMach012 Generator model from OpenDSS.

Based on the following files from the official OpenDSS source code:
- Source/PCElements/IndMach012.pas
- Source/IndMach012a/IndMach012Model.pas

This Python version was written by Paulo Meira. 
Original code by EPRI, licensed under the 3-clause BSD. See OPENDSS_LICENSE.

This sample code doesn't interact with the main OpenDSS interface directly,
it only uses the user-model interface. Thus, it is compatible with the official OpenDSS
distribution as well as DSS-Python. Note that OpenDSS version 7 has a bug on 64-bit 
system and user-models most likely won't run via COM.

Recent version of OpenDSS 8 also present a bug when handling the edition of 
user-model parameters after the creation of the generator. You can, of course,
edit the data in Python if you desire.

from dss.UserModels import GenUserModel
from dss.enums import SolveModes
import numpy as np

# This is the user-model we'll use as a base
Base = GenUserModel.Base

# Symmetrical component transformation matrices
a = np.exp(1j * 2 * np.pi/3)
aa = np.exp(1j * 4 * np.pi/3)

Ap2s = np.array([
    [1, 1, 1],
    [1, a, aa],
    [1, aa, a]
]) / 3.0

As2p = np.array([
    [1, 1, 1],
    [1, aa, a],
    [1, a, aa]

@GenUserModel.register # The class needs to be registered
class PyIndMach012(Base):
    A Python User-Model implementation of the IndMach012 Generator model for OpenDSS 
    def __init__(self, gen, dyn, callbacks):
        Initialize the model object instance
        Note: OpenDSS calls `Init` from the UserModel DLL later,
              which calls our `init_state_vars`

        Base.__init__(self, gen, dyn, callbacks)

        # You can list the DSS model inputs and default values like this.
        # The UserModel wrapper will create attributes in the instance with
        # the default values and update them later when the UserModel
        # `Edit` function is called.
            ('H', 0.02),
            ('D', 0.02),
            ('puRs', 0.0053),
            ('puXs', 0.106),
            ('puRr', 0.007),
            ('puXr', 0.12),
            ('puXm', 4.0),
            ('slip', 0.007),
            ('MaxSlip', 0.1),
            ('slipOption', 'variable')

        # The outputs can be any variable or Python property, i.e. it can be an 
        # input, state variable, property, etc., as long as it is available in
        # the model class
            'Slip', # The current slip (`slip` is the DSS input param)

            # There don't need to be in the output (they're constant) but are listed 
            # in IndMach012.pas -- most likely to debug

            # complex variables like these are exported as their absolute values

            # Some properties to mimic the Pascal version

        # These are the state variables. DSS-Python will automatically
        # setup auxiliary variables such as dE1_dt, dE1_dtn, E1n used in
        # the solution process
            'E1', 'E2'

        # For some advanced usage, we need some CFFI code.
        # We plan to add a simple wrapper to the callback interface.
        # While writing this, I noticed that there were some changes in 
        # OpenDSS version 8 that introduce an extra ActorID parameter in 
        # many of the callback functions. This means that we cannot write
        # a user-model that is compatible with both versions. 
        # An issue ticket was created to track this at:
        # self.char_buffer ='char[1024]')
        # self.callbacks.GetActiveElementName(self.char_buffer, 1024)
        # self.element_name = self.ffi.string(self.char_buffer)#.decode('ascii')

        # This one is used for the Power property, left as an example
        # self.double_buffer ='double[2]')

        # Update other variables that depend on the input parameters

    def init_state_vars(self, Vabc, Iabc):
        Initialize state variables (dynamics mode), equivalent to

        V012 =, Vabc)
        I012 =, Iabc)

        # The following is done in TIndMach012Obj.InitModel:
        # Compute Voltage behind transient reactance and set derivatives to zero
        self.E1 = V012[1] - I012[1] * self.Zsp
        self.dE1dt = 0
        self.E2 = V012[2] - I012[2] * self.Zsp
        self.dE2dt = 0

        # Copy the current state to the previous state

        # Initial rotor speed
        self.gen.Speed = -self.S1 * self.gen.w0
        self.gen.dSpeed = 0
        self.gen.Theta = np.angle(self.E1) # overwrite Theta
        self.gen.dTheta = 0

    def integrate(self):
        Equivalent to TIndMach012Obj.Integrate
        if self.dyn.IterationFlag == 0:
            # First iteration of new time step, copy the previous state
            # to be used in the integration process

        # Some copies to reduce `self.` spam
        gen = self.gen
        w0 = gen.w0
        S1, S2 = self.S1, self.S2
        E1, E2 = self.E1, self.E2
        Is1, Is2 = self.Is1, self.Is2
        T0p = self.T0p
        Xopen, Xp = self.Xopen, self.Xp

        # Derivative of E
        self.dE1_dt = (1j * -w0 * S1 * E1) - ((E1 - 1j * (Xopen - Xp) * Is1) / T0p)
        self.dE2_dt = (1j * -w0 * S2 * E2) - ((E2 - 1j * (Xopen - Xp) * Is2) / T0p)

        # Trapezoidal Integration

    def update(self):
        Propagate changes from the input parameters to the model.
        Equivalent to part of TIndMach012Obj.RecalcElementData

        gen = self.gen


        # make generator speed agree
        gen.Speed = -self.S1 * self.gen.w0
        self.gen.dSpeed = 0.0

        self.fixed_slip = (self.slipOption) and (self.slipOption[0].upper() == 'F')
        self.first_iteration = True

        ZBase = 1000.0 * (gen.kVGeneratorBase**2 / gen.kVArating)

        Rs = self.puRs * ZBase
        Xs = self.puXs * ZBase
        Rr = self.puRr * ZBase
        Xr = self.puXr * ZBase
        Xm = self.puXm * ZBase

        self.Zs = complex(Rs, Xs)
        self.Zm = complex(0, Xm)
        self.Zr = complex(Rr, Xr)

        self.Xopen = Xs + Xm
        self.Xp  = Xs + (Xr * Xm) / (Xr + Xm)
        self.Zsp = complex(Rs, self.Xp)
        self.T0p = (Xr + Xm) / (gen.w0 * Rr)
        # self.Zrsc = self.Zr + (self.Zs * self.Zm) / (self.Zs + selfg.Zm)

        # Init dSdP based on rated slip and rated voltage
        self.V1 = complex(gen.kVGeneratorBase * 1000.0/np.sqrt(3))
        if self.S1 != 0:
            self.Is1, self.Ir1 = self._pfmodel_current(self.V1, self.S1)

        self.dSdP = self.S1/(self.V1 * np.conjugate(self.Is1)).real

        self.Is1 = complex(0)
        self.V1 = complex(0)
        self.Is2 = complex(0)
        self.V2 = complex(0)

    def calc(self, Vabc, Iabc):
        Calculate the new model state. Vabc is used as an
        input, while Iabc is the ouput used in OpenDSS.

        # The next version of DSS-Python should have an option to 
        # provide the values in 012 space to simplify the model code
        V012 =, Vabc)
        I012 =, Iabc)

        if self.dyn.SolutionMode == SolveModes.Dynamic:
            self.calc_dynamic(V012, I012)
            self.calc_power_flow(V012, I012)

        Iabc[:] = iabc =, I012)

        # We can get the current total power here, or we can use 
        # the power property below 
        self.Power = sum(np.asarray(Vabc) * iabc.conj())

    # @property
    # def Power(self):
        # '''
        # This is an example of callback usage, returning the power of the 
        # element. Note that we don't really need this here since we can 
        # calculate the power in the calc function.
        # '''
        # cmd = b'select %s' % (self.element_name)
        # self.callbacks.DoDSSCommand(cmd, len(cmd) + 1)
        # self.callbacks.GetActiveElementPower(1, self.double_buffer)
        # return complex(self.double_buffer[0], self.double_buffer[1])

    def calc_dynamic(self, V012, I012):
        '''Equivalent to TIndMach012Obj.CalcDynamic'''

        # In dynamics mode, slip is allowed to vary

        # Copy some values to local variables
        V1, V2 = self.V1, self.V2 = V012[1], V012[2]
        E1, E2 = self.E1, self.E2
        Zsp, Zm = self.Zsp, self.Zm

        # Gets slip from shaft speed
        self._set_local_slip(-self.gen.Speed / self.gen.w0)

        # The stator and rotor currents from the Pascal code are 
        # computed in TIndMach012Obj.Get_DynamicModelCurrent

        # Stator current
        self.Is1 = (V1 - E1) / self.Zsp
        self.Is2 = (V2 - E2) / self.Zsp

        # Rotor current
        self.Ir1 = self.Is1 - (V1 - self.Is1 * Zsp) / Zm
        self.Ir2 = self.Is2 - (V2 - self.Is2 * Zsp) / Zm

        I012[:] = complex(0.0, 0.0), self.Is1, self.Is2

    def calc_power_flow(self, V012, I012):
        '''Equivalent to TIndMach012Obj.CalcPFlow'''

        self.V1, self.V2 = V012[1], V012[2]

        if self.first_iteration:
            # Initialize Is1
            self.Is1, self.Ir1 = self._pfmodel_current(self.V1, self.S1)

        # If fixed slip option set, then use the value set by the user
        if not self.fixed_slip:
            P_Error = self.gen.PNominalPerPhase - (self.V1 * self.Is1.conjugate()).real

            # make new guess at slip
            self._set_local_slip(self.S1 + self.dSdP * P_Error)

        self.Is1, self.Ir1 = self._pfmodel_current(self.V1, self.S1)
        self.Is2, self.Ir2 = self._pfmodel_current(self.V2, self.S2)

        I012[:] = complex(0.0, 0.0), self.Is1, self.Is2

    def _pfmodel_current(self, V, s, show=False):
        '''Equivalent to TIndMach012Obj.Get_PFlowModelCurrent'''

        if s != 0.0:
            RL = self.Zr.real * (1 - s) / s
            RL = self.Zr.real * 1.0e6

        Zrotor = RL + self.Zr
        Zmotor = self.Zs + (Zrotor * self.Zm) / (Zrotor + self.Zm)
        Istator = V / Zmotor
        Irotor = Istator - (V - self.Zs * Istator) / self.Zm

        return Istator, Irotor

    def _set_local_slip(self, value):
        '''Equivalent to TIndMach012Obj.set_Localslip'''

        self.S1 = value
        if self.dyn.SolutionMode != SolveModes.Dynamic:
            # Put limits on the slip unless dynamics
            if abs(self.S1) > self.MaxSlip:
                self.S1 = np.sign(self.S1) * self.MaxSlip

        self.S2 = 2 - self.S1

    # The following are properties to emulate the model outputs from 
    # the Pascal version of built-in IndMach012 

    def E1_pu(self):
        return np.sqrt(3) * abs(self.E1) / (1000 * self.gen.kVGeneratorBase)

    def Slip(self):
        return self.S1

    def RotorLosses(self):
        Ir1, Ir2, Zr = self.Ir1, self.Ir2, self.Zr
        return 3 * (Ir1.real**2 + Ir1.imag**2 + Ir2.real**2 + Ir2.imag**2) * Zr.real

    def StatorLosses(self):
        Is1, Is2, Zs = self.Is1, self.Is2, self.Zs
        return 3 * (Is1.real**2 + Is1.imag**2 + Is2.real**2 + Is2.imag**2) * Zs.real

    def PowerFactor(self):
        power = self.Power
        return np.sign(power.imag) * power.real / abs(power)

    def Efficiency_pct(self):
        power = self.Power
        return np.clip((1 - (self.StatorLosses + self.RotorLosses) / power.real) * 100, 0, 100)

    def ShaftPower_hp(self):
        Ir1, Ir2, Zr, S1, S2 = self.Ir1, self.Ir2, self.Zr, self.S1, self.S2
        return (3.0/746) * (abs(Ir1)**2 * (1 - S1) / S1 + abs(Ir2)**2 * (1 - S2)/S2) * Zr.real

OpenDSS setup

For this example, we can use either COM or DSS-Python (DSS C-API).

original_dir = os.getcwd() # same the original working directory since the COM module messes with it
USE_COM = False # toggle this value to run with DSS C-API or COM
    from dss import patch_dss_com
    import comtypes.client
    DSS = patch_dss_com(comtypes.client.CreateObject('OpenDSSengine.DSS'))
    DSS.DataPath = original_dir
    from dss import DSS

'DSS C-API Library version DEV revision UNKNOWN based on OpenDSS SVN UNKNOWN [FPC 3.2.2] (64-bit build) MVMULT INCREMENTAL_Y CONTEXT_API PM 20230330230049; License Status: Open '
Text = DSS.Text
Monitors = DSS.ActiveCircuit.Monitors

Using the model

To use a Python model for generators: - the model class needs to be registered in advance - create a generator with model=6 - pass a usermodel="{dll_path}" as in the following DSS command in the run function - pass a "pymodel=MODELNAME" parameter in the userdata property, where MODELNAME is the name of the model class in Python

def run(pymodel):
    Text.Command = 'redirect "master.dss"'

    if pymodel:
        # This uses our custom user-model in Python
        Text.Command = 'New "Generator.Motor1" bus1=Bg2 kW=1200 conn=delta kVA=1500.000 H=6 model=6 kv=0.48 usermodel="{dll_path}" userdata=(pymodel=PyIndMach012 purs=0.048 puxs=0.075 purr=0.018 puxr=0.12 puxm=3.8 slip=0.02 SlipOption=variableslip)'.format(
        Text.Command = 'New "Monitor.mfr2" element=Generator.Motor1 terminal=1 mode=3'
        # This uses the built-in model for comparison
        Text.Command = 'New "IndMach012.Motor1" bus1=Bg2 kW=1200 conn=delta kVA=1500.000 H=6 purs=0.048 puxs=0.075 purr=0.018 puxr=0.12 puxm=3.8 slip=0.02 SlipOption=variableslip kv=0.48'
        Text.Command = 'New "Monitor.mfr2" element=IndMach012.Motor1 terminal=1 mode=3'

    # This will run a power-flow solution
    Text.Command = 'Solve'

    # This will toggle to the dynamics mode
    Text.Command = 'Set mode=dynamics number=1 h=0.000166667'

    # And finally run 5000 steps for the dynamic simulation
    Text.Command = f'Solve number=5000'
# There are the channels from the Pascal/built-in IndMach012
channels_pas = ('Frequency', 'Theta (deg)', 'E1', 'dSpeed (deg/sec)', 'dTheta (deg)', 'Slip', 'Is1', 'Is2', 'Ir1', 'Ir2', 'Stator Losses', 'Rotor Losses', 'Shaft Power (hp)', 'Power Factor', 'Efficiency (%)')

# There are the channels from the Python module -- we define part of these and part come from the generator model itself
channels_py = ('Frequency', 'Theta (Deg)', 'E1_pu', 'dSpeed (Deg/sec)', 'dTheta (Deg)', 'Slip', 'Is1', 'Is2', 'Ir1', 'Ir2', 'StatorLosses', 'RotorLosses', 'ShaftPower_hp', 'PowerFactor', 'Efficiency_pct')

Running and saving the outputs

Let’s run the Pascal/built-in version of IndMach012 and our custom Python version for comparison:

Monitors.Name = 'mfr2'
header = [x.strip() for x in Monitors.Header]
outputs_pas = {channel: Monitors.Channel(header.index(channel) + 1) for channel in channels_pas}

Monitors.Name = 'mfr2'
header = [x.strip() for x in Monitors.Header]
outputs_py = {channel: Monitors.Channel(header.index(channel) + 1) for channel in channels_py}

time = np.arange(1, 5000 + 1) * 0.000166667
offset = int(0.1 / 0.000166667)

Plotting the various output channels

The example circuit applies a fault at 0.3 s, isolating the machine at 0.4s (check master.dss for more details).

As we can see from the figures below, the outputs match very closely. After the induction machine is isolated, the efficiency and power factor values can misbehave as the power goes to zero, seem especially in the Pascal version.

for ch_pas, ch_py in zip(channels_pas, channels_py):
    plt.plot(time, outputs_pas[ch_pas], label='IndMach012', lw=3)
    plt.plot(time, outputs_py[ch_py], label='PyIndMach012', ls='--', lw=2)
    plt.axvline(0.3, linestyle=':', color='k', alpha=0.5, label='Fault occurs')
    plt.axvline(0.4, linestyle='--', color='r', alpha=0.5, label='Relays operate')
    plt.xlabel('Time (s)')

    if ch_pas == 'Efficiency (%)':
        # Limit efficiency to 0-100
        plt.ylim(0, 100)

    plt.xlim(0, time[-1])