Modify Smooth Approximation of Limits

[lukelowry]

Description

Change the smooth limit approximation used in the IEEET1 Exciter and TGOV1 Governor so that it is twice differentiable. Additionally, the exciter saturation equation now has a smooth approximation. This is for Enzyme compatibility.

Proposed changes

The logistic function is now used instead of the approximation of the logistic function

$$\dfrac{\alpha x}{1+|\alpha x|} \to \dfrac{1}{1+\exp(-\alpha x)}$$

Checklist

• All tests pass.
• Code compiles cleanly with flags -Wall -Wpedantic -Wconversion -Wextra.
• The new code follows GridKit™ style guidelines.
• There are unit tests for the new code.
• The new code is documented.
• The feature branch is rebased with respect to the target branch.
• N/A: I have updated CHANGELOG.md to reflect the changes in this PR. If this is a minor PR that is part of a larger fix already included in the file, state so.

Further Comments

Needed for Enzyme to work on larger cases

[nkoukpaizan]

@lukelowry Thanks for this. I likely didn't make myself clear in the meeting last week. The challenging part for Enzyme is not the sigmoid, but the saturation: the following branching in the exciter model.

GridKit/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp

Line 346 in 6202f80

if (efdp > SA_)

[lukelowry]

@lukelowry Thanks for this. I likely didn't make myself clear in the meeting last week. The challenging part for Enzyme is not the sigmoid, but the saturation: the following branching in the exciter model.

Ah, I see now. Fixed in the recent commit. Please advise if this is sufficient for Enzyme.

[nkoukpaizan]

Looks like the new sigmoid is more challenging to differentiate (sharper approximation?) than the previous one.

One of the derivatives is nan for Tgov1 with both Enzyme and the DependencyTracking chain rule approaches

46: 0th residual: 0 dependencies: [ (0, -0.133333) (1, 0.666667) ]
46: 1th residual: 0 dependencies: [ (1, nan) (3, -30) ]
46: 2th residual: 0 dependencies: [ (0, 0.133333) (1, 0.333333) (2, -1) (3, 0) ]
46: Evaluate Jacobian for Tgov1...
46: Jacobian evaluation is experimental!
46: Sparse COO Matrix: Model Jacobian
46: (x , y, value)
46: (0, 0, -0.133333)
46: (0, 1, 0.666667)
46: (1, 1, nan)
46: (1, 3, -30)
46: (2, 0, 0.133333)
46: (2, 1, 0.333333)
46: (2, 2, -1)
46: (2, 3, -0)

The previous result was:

46: 0th residual: 0 dependencies: [ (0, -0.133333) (1, 0.666667) ]
46: 1th residual: 0 dependencies: [ (1, -1.49991) (3, -29.9981) ]
46: 2th residual: 0 dependencies: [ (0, 0.133333) (1, 0.333333) (2, -1) (3, 0) ]
46: Evaluate Jacobian for Tgov1...
46: Jacobian evaluation is experimental!
46: Sparse COO Matrix: Model Jacobian
46: (x , y, value)
46: (0, 0, -0.133333)
46: (0, 1, 0.666667)
46: (1, 1, -1.49991)
46: (1, 3, -29.9981)
46: (2, 0, 0.133333)
46: (2, 1, 0.333333)
46: (2, 2, -1)
46: (2, 3, -0)

[pelesh]

Looks like the new sigmoid is more challenging to differentiate (sharper approximation?) than the previous one.

Just take a look at parameters -- the code is using $\dfrac{1}{1+\exp(-\alpha x)}$ with $\alpha = 1000$. For numerical purposes, this is a discrete switch.

[pelesh]

@lukelowry, what is the lowest value of $\alpha$ ( in code it is mu), which you could use in this approximation without jeopardizing the accuracy?

[lukelowry]

@pelesh @nkoukpaizan Hard to say what the threshold is, but since the limits are on the order of $10^0$, $\alpha\approx 50$ is the lowest I would set it.

We can change it back to the absolute value form, if needed.

[pelesh]

@pelesh @nkoukpaizan Hard to say what the threshold is, but since the limits are on the order of 10 0 , α ≈ 50 is the lowest I would set it.

We can change it back to the absolute value form, if needed.

I think exponential form is better, one just need to be mindful of the exponential growth rate. $\alpha$ is not really an arbitrary parameter but $\alpha \sim 1/\tau$ where $\tau$ is the smallest time scale you want to resolve in your simulations.

[nkoukpaizan]

From the derivative's perspective, the function is well-behaved up to $\alpha \sim 700$ in the exponential form. All other simulation examples pass, but not sure if the switch is truly being put to the test in these examples.

[nkoukpaizan]

I have confirmed that removing the branching gets me past the Enzyme build error I was seeing. I'm seeing another error that is unrelated to the model implementation, and I can work around that.

For this PR, I would just suggest setting a lower $\alpha$ value (e.g., $4 \times 60=240$, matching $1/\Delta t$ in examples) with the exponential form of the logistic function.

[nkoukpaizan]

Since I made the latest changes, it would be good to have another pair of eyes (@pelesh @abirchfield ) on this before merging.

[pelesh]

For some reason exciter tests fail for me:

The following tests FAILED:
         25 - GenrouTest1_Ieeet1_Json (Subprocess aborted)
         26 - GenrouTest1_Ieeet1_Json_no_arg (Subprocess aborted)
Errors while running CTest
Output from these tests are in: /Users/55y/src/gridkit/build-gcc/Testing/Temporary/LastTest.log

[pelesh]

I am not sure how this passed CI pipeline:

26: Test timeout computed to be: 1500
26: Example: TwoBusIeeet1Json
26: Input file: "/Users/55y/src/gridkit/build-gcc/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.json"
26: [ERROR] Genrou: pmech signal attached with no linked governor
26: [ERROR] Component errors: 1
26: terminate called after throwing an instance of 'std::runtime_error'
26:   what():  SystemModel allocation failed
2/2 Test #26: GenrouTest1_Ieeet1_Json_no_arg ...Subprocess aborted***Exception:   0.03 sec

[pelesh]

@PhilipFackler, I believe the check here may be overzealous:

GridKit/GridKit/Model/PhasorDynamics/SystemModel.hpp

Line 332 in 8c92cc2

if (errorCount > 0)

If I comment out the exception, everything works correctly.

I don't think that at this point we can check if a governor is connected to the generator because governor may have not been allocated/connected at this point. Also, we may need to discriminate here between sensors and actuators.

CC @nkoukpaizan