Question: behavior of sigmoid and Indicators

[nkoukpaizan]

A couple of thoughts on the sigmoid and indicator functions in IEEET1 and Tgov1.

1. In Modify Smooth Approximation of Limits #313, I moved the definition of these math functions to a common file, such that all models will use the same functions. However, I am finding out in testing PhasorDynamics simulations with sparse Jacobians #312 that the values of $\mu$ for finite numerical derivatives is not the same for IEEET1 and Tgov1. I need a lower value of $\mu$ for the IEEET1 example after the fault ($\sim 12$ versus $240$). To me, this points to the need to normalize the variables in the exponential form $\left( \dfrac{1}{1+e^{-\alpha x}} \to \dfrac{1}{1+e^{-\alpha x/x_{ref}}} \right)$ to ensure the derivative is finite. Otherwise, it will always be possible to pass values of $x$ that will break the derivative, even if we set $\mu$ to lower values. The alternative would be to revert to the absolute value form, whose derivative is less sensitive.
2. The behavior of the indicator is significantly different when $f \simeq 0$ versus $\lvert f\rvert >> 0$ and I am concerned this will affect the verification with other tools. In particular below is a plot with $\sigma(f)\rvert_{0} =0.5$, where the indicator values are 0.5 outside of the target range. I've included $\sigma(f)\rvert_{-0.25} = 0$ and $\sigma(f)\rvert_{0.25} = 1$ results as well for comparison. If I read the equations correctly, the desired value outside of the target range is 0.

$\sigma(f)\rvert_{-0.25} = 0$

$\sigma(f)\rvert_{0} =0.5$

$\sigma(f)\rvert_{0.25} = 1$

Below is the code I used to plot, with the parameter choices

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import interact
# Enable the matplotlib widget backend in a Jupyter environment
%matplotlib widget 

plt.rc('text', usetex=False)
plt.rc('font', family='serif')
plt.rc('figure', figsize=(12,5))
opacity   = 0.9
linewidthVal = 2.0
fontsizeVal = 16
colors = ['k', 'b', 'r', 'c', 'y', 'm', 'g']

def sigmoid_exp(x):
    MU = 240.0
    return 1.0 / (1.0 + np.exp(-MU * x))

def sigmoid_abs(x):
    MU = 40000
    return (0.5 * MU * x) / (1.0 + np.abs(MU * x)) + 0.5

def sigmoid(x):
    return sigmoid_exp(x)

def indicator_low(limit_min, x, f):
    return sigmoid(limit_min - x) * sigmoid(-f)

def indicator_high(limit_max, x, f):
    return sigmoid(x - limit_max) * sigmoid(f)

def indicator(limit_min, limit_max, x, f):
    return (1.0 - indicator_low(limit_min, x, f)) * (1.0 - indicator_high(limit_max, x, f))

def make_plot(func):
    x_min = -0.5
    x_max = 0.5
    x = np.linspace(2*x_min, 2*x_max, 1000)
    
    plt.clf()
    plt.plot(x, sigmoid_exp(x), '-', linewidth=linewidthVal, color=colors[0], label='sigmoid_exp (240)')
    plt.plot(x, indicator_low(x_min, x, func), '-', linewidth=linewidthVal, color=colors[1], label='indicator_low')
    plt.plot(x, indicator_high(x_max, x, func), '-', linewidth=linewidthVal, color=colors[2], label='indicator_high')
    plt.plot(x, indicator(x_min, x_max, x, func), '-', linewidth=linewidthVal, color=colors[3], label='indicator')
    plt.xlabel('x', fontsize=fontsizeVal*1.2, style='italic')
    plt.ylabel('y', fontsize=fontsizeVal*1.2, style='italic')
    plt.legend(fontsize=fontsizeVal)
    plt.xticks(fontsize=fontsizeVal)
    plt.yticks(fontsize=fontsizeVal)
    plt.show()

interact(make_plot, func=(-1.0, 1.0, 0.25));

Thoughts?
CC @pelesh @abirchfield @lukelowry

[pelesh]

Yes, in these functions $x$ should come in already normalized, if we are to have a system-wide $\alpha$. We should probably add that to the modeling guidelines. We could, for example, require that a model developer provide $x_{ref}$ for each variable.

Note that most of the phasor domain variables and parameters are already normalized in terms of "power units".

[nkoukpaizan]

For the second point I raised regarding the behavior for $f \simeq 0$, we could make use of a normalized approximation to the Dirac delta function $\delta(f) = 4 \sigma(f) \left(1.0 - \sigma(f) \right)$ using the derivative of the sigmoid to be explicit about the $f=0$ case.

The indicator would be of the form $\phi = \phi_0 \delta(f) + \left(\phi_L \sigma(-f) + \phi_U \sigma(f) \right) \left(1- \delta(f) \right)$,
where
$\phi_L = \sigma(x - x_{min}) = 1 - \sigma(x_{min} - x) $,
$\phi_U = \sigma(x_{max} - x) = 1 - \sigma(x - x_{max} )$
$\phi_0 = \phi_L + \phi_U - 1$

This would result in the following behavior when $f = 0$:

This does make the indicator worse for small non-zero values of $f$. Basically interpolating in [0, 1] versus [0.5, 1] converges at higher absolute values of $f$. e.g., $f = -0.01$:

There is no change for $\lvert f \rvert >> 0$

Below is the code I used to plot, with the parameter choices

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import interact
# Enable the matplotlib widget backend in a Jupyter environment
%matplotlib widget 

plt.rc('text', usetex=False)
plt.rc('font', family='serif')
plt.rc('figure', figsize=(12,5))
opacity   = 0.9
linewidthVal = 2.0
fontsizeVal = 16
colors = ['k', 'b', 'r', 'c', 'y', 'm', 'g']

def sigmoid_exp(x):
    MU = 300.0
    return 1.0 / (1.0 + np.exp(-MU * x))

def sigmoid_abs(x):
    MU = 4000.0
    return (0.5 * MU * x) / (1.0 + np.abs(MU * x)) + 0.5

def sigmoid(x):
    return sigmoid_exp(x)

def dsigmoid(x): # normalized for unitary peak value
    return 4 * sigmoid(x) * (1.0 - sigmoid(x))

def indicator_low(limit_min, x, f):
    return sigmoid(limit_min - x) * sigmoid(-f)

def indicator_high(limit_max, x, f):
    return sigmoid(x - limit_max) * sigmoid(f)

def indicator(limit_min, limit_max, x, f):
    return (1.0 - indicator_low(limit_min, x, f)) * (1.0 - indicator_high(limit_max, x, f))

def my_indicator2(limit_min, limit_max, x, f):
    df = dsigmoid(f)
    low = sigmoid(x - limit_min)
    high = sigmoid(limit_max - x)
    zero = low + high - 1
    return zero * df + (low * sigmoid(-f) + high * sigmoid(f)) * (1.0 - df)

def make_plot(func):
    x_min = -0.5
    x_max = 0.5
    x = np.linspace(2*x_min, 2*x_max, 1000)
    
    plt.clf()
    plt.plot(x, sigmoid_exp(x), '-', linewidth=linewidthVal, color=colors[0], label='sigmoid_exp (240)')
    #plt.plot(x, sigmoid_abs(x), '--', linewidth=linewidthVal, color=colors[0], label='sigmoid_abs (40000)')
    plt.plot(x, dsigmoid(x), '-.', linewidth=linewidthVal, color=colors[0], label='dsigmoid')
    #plt.plot(x, indicator_low(x_min, x, func), '-', linewidth=linewidthVal, color=colors[1], label='indicator_low')
    #plt.plot(x, indicator_high(x_max, x, func), '-', linewidth=linewidthVal, color=colors[2], label='indicator_high')
    plt.plot(x, indicator(x_min, x_max, x, func), '-', linewidth=linewidthVal, color=colors[3], label='old indicator')
    plt.plot(x, my_indicator2(x_min, x_max, x, func), '-', linewidth=linewidthVal, color=colors[4], label='new indicator')
    plt.xlabel('x', fontsize=fontsizeVal*1.2, style='italic')
    plt.ylabel('y', fontsize=fontsizeVal*1.2, style='italic')
    plt.legend(fontsize=fontsizeVal)
    plt.xticks(fontsize=fontsizeVal)
    plt.yticks(fontsize=fontsizeVal)
    plt.show()

interact(make_plot, func=(-0.05, 0.05, 0.01));
#make_plot(0.5)