Phasor Dynamics examples do not re-solve for consistent initial conditions after discrete events, such as bus faults

[alexander-novo]

Overview

As an example, consider the PhasorDynamics/TenGen example, which uses bus faults as discrete events during the simulation. After (de)activating the bus fault, the simulation is re-initialized:

GridKit/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp

Lines 172 to 175 in b191e66

	fault.setStatus(1);
	ida.initializeSimulation(1.0, false);
	nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
	ida.runSimulation(1.1, nout, output_cb);

But a value of false is passed to to initializeSimulation for the parameter findConsistent, which is used to determine if IDACalcIC() should be called:

GridKit/src/Solver/Dynamic/Ida.cpp

Lines 270 to 282 in b191e66

	if (findConsistent)
	{
	int initType = IDA_Y_INIT;
	if (tag_)
	initType = IDA_YA_YDP_INIT;
	retval = IDACalcIC(solver_, initType, 0.1);
	checkOutput(retval, "IDACalcIC");
	copyVec(yy_, model_->y());
	copyVec(yp_, model_->yp());
	}

The documentation for IDACalcIC() says

IDACalcIC() calculates corrected initial conditions for the DAE system for certain index-one problems including a class of systems of semi-implicit form (see §5.2.2 and [26]). It uses a Newton iteration combined with a linesearch algorithm. Calling IDACalcIC() is optional. It is only necessary when the initial conditions do not satisfy the given system. Thus if y0 and yp0 are known to satisfy $F(t_0, y_0, \dot{y}_0) = 0$, then a call to IDACalcIC() is generally not necessary.

In general, we don't know that the condition $F(t_0, y_0, \dot{y}_0) = 0$ is satisfied, and by experimentation we can verify that it isn't in this particular example. This indicates that a call to IDACalcIC() should still be necessary.

The documentation also says

A call to the function IDACalcIC() must be preceded by successful calls to IDACreate() and IDAInit() (or IDAReInit()), and by a successful call to the linear system solver specification function. The call to IDACalcIC() should precede the call(s) to IDASolve() for the given problem.

This seems to imply to me that a call to IDAReInit() (which is being done in initializeSimulation()) is not enough to verify that initial conditions are consistent and that we must still call IDACalcIC() (@Steven-Roberts can you confirm?)

Background

I've known about this for some time - in our Matlab models which are being compared to solutions generated by GridKit, the values after discrete events would always differ drastically, and step size controllers of certain Matlab integrators were not able to meet requested tolerances when we would not re-solve for consistent initial conditions (not all integrators are as tolerant as BDF turns out haha). When I discussed it with @reid-g, he mentioned that this was expected behavior - power systems engineers would not want algebraics to change discretely when re-solving for consistent initial conditions. But when I let @abirchfield and Carol know in our meeting about fixed vs adaptive integrators, they were surprised to learn about this behavior.

Fixes

If this behavior is intentional, it (and the reasoning behind it) needs to be documented somewhere since multiple people now have been surprised by it. Also maybe we should have a discussion with @abirchfield and Carol verifying that's what we want.

Otherwise, it would be best to fix these examples by changing the value false to true. Also, it would be interesting to see if this API could be changed so that true is the default - it doesn't happen often in our simulations, I think, that the simulation needs to be re-initialized and most of those times the initial conditions should be made consistent. The extra (potentially unneeded) computation should be overshadowed by the rest of the simulation.

[Steven-Roberts]

This seems to imply to me that a call to IDAReInit() (which is being done in initializeSimulation()) is not enough to verify that initial conditions are consistent and that we must still call IDACalcIC() (@Steven-Roberts can you confirm?)

Correct, IDAReInit neither checks for nor enforces consistency. If the inconsistency is on the order of the tolerances, this is usually benign. For larger inconsistencies, I would recommend calling IDACalcIC.

I agree it would be safest to always solve for consistent initial conditions after re-initialization. It's also very important to ensure initType is set correctly. After a fault, I think IDA_YA_YDP_INIT makes the most sense.

[abirchfield]

I am OK with adding the IDACalcIC call. However, we would need to investigate if this change caused our currently validated examples to change significantly.

[pelesh]

I think the proposed feature is important and should be added to GridKit's configuration of the numerical solver.

For the models we currently have in GridKit, this is not an urgent issue, though. As far as I understand SUNDIALS code, if you call IDACalcIC, you will be using a linesearch Newton method. If you use reinitialization only, you would be using "lightweight" Newton solver that is embedded in the IDA numerical integrator. The difference between the two is that linesearch Newton has larger basin of convergence than the lightweight variant.

For time-triggered events, IDACalcIC is typically not needed; the lightweight solver will take care of it in most cases. On the other hand, for state-triggered events one will likely need to use IDACalcIC. Once we develop components that induce state-triggered events, we would need to add an option to use search for consistent set of initial conditions upon the event.

[alexander-novo]

@pelesh Sorry for my late reply to this, I misunderstood what you meant originally, but I now think I understand what you mean.

The newton's method used in IDACalcIC should only solve for the algebraic variables of the state, but the newton's method used inside the method should solve for the entire state, including differentials. As well, the newton's method inside the method should be used to advance the state - i.e. it will be used to calculate what the state should look like some time in the future - I don't think it will be used to calculate anything at the current time (i.e. initial conditions).

The convergence of the numerical scheme relies on the consistency of the initial conditions - error that exists in the initial conditions may propagate to the final solution.

Please let me know if my understanding of what's going on in IDA is consistent @Steven-Roberts

[Steven-Roberts]

@alexander-novo yes, that's correct.

If there is any chance that an event would change the simulation state (or derivative), I'd recommend calling IDACalcIC after IDAReInit. I'm worried that if we try to optimize away calls to IDACalcIC we risk some subtle bugs when the state changes more than we anticipated. If IDACalcIC it happens to be unnecessary for a particular event, then that will cost a couple extra function evaluation and one linear solve. Depending on the frequency of events, this could be negligible.

[alexander-novo]

I also think the IDA_Y_INIT mode (which is the default in Ida.cpp) is never what we want when solving for consistent initial conditions. We should always use IDA_YA_YDP_INIT mode, and unavailability of a tag vector should be an error. The use of IDA_Y_INIT combined with the implementation of COO_matrix::axpy is the reason for the segmentation fault I mentioned in today's meeting @pelesh (and even fixing this by fixing COO_matrix::axpy will produce numerical errors).