Initial conditions solving.

[Federerer]

So, during the work on dynamic parameters, I've noticed that initial condition solver stopped working. In my current solution I'm adding an unknown for both mapped value of the potentiometer (with taper function applied), top resistance and bottom resistance, as I've noticed that it helps with not getting insanely large solutions, as now all the resistance calculation ends up in the linear part of the solution. This unfortunately broke the steady state solver because it sets initial state of the unknowns to zero:

List<Arrow> part = i.Equations.Select(j => Equal.New(j, 0)).NSolve(i.Unknowns.Select(j => Arrow.New(j, 0)));

this is fine for voltages, but not for resistances 😁 So, I tried to solve the linear part first, but I've noticed that for some circuits back substitution takes very long time and produces huge expressions, and I don't need to back substitute the whole system except the linear part. I'm aware that I can provide a list of unknowns to BackSubstitute method, but this is sort of a chicken-egg problem, because I don't know the solvable unknowns upfront. So I came up with this ugly thing, that does only a row reduction and creates a new system just for linear variables:

// Find steady state solution for initial conditions.
List<Arrow> initial = InitialConditions.ToList();
Log.WriteLine(MessageType.Info, "Performing steady state analysis...");

var defaultArguments = Analysis.Parameters.Select(i => Arrow.New(i.Expression, i.Value)).ToList();

SystemOfEquations dc = new SystemOfEquations(equations
    // Substitute default arguments for the parameters.
    .Evaluate(defaultArguments)
    // Derivatives, t, and T are zero in the steady state.
    .Substitute(dy_dt.Select(i => Arrow.New(i, 0)).Append(Arrow.New(t, 0), Arrow.New(T, 0)))
    // Use the initial conditions from analysis.
    .Substitute(Analysis.InitialConditions)
    // Evaluate variables at t=0.
    .OfType<Equal>(), unknowns.Select(j => j.Substitute(t, 0)));

// Solve partitions independently.
foreach (SystemOfEquations i in dc.Partition())
{
    try
    {
        i.RowReduce();
        //i.BackSubstitute(); // insanely slow for some circuits
        var lin = i.Solve();

        var linearSystem = new SystemOfEquations(lin.Select(l => Equal.New(l.Right, l.Left)), lin.Select(l => l.Left));
        linearSystem.RowReduce();
        linearSystem.BackSubstitute();
        var x = linearSystem.Solve();
        initial.AddRange(x);

        var other = i.Evaluate(x);
        var init = i.Unknowns.Select(j => Arrow.New(j, 0)).ToArray();
        List<Arrow> part = other.Select(j => Equal.New(j, 0)).NSolve(init);
        initial.AddRange(part);
        LogExpressions(Log, MessageType.Verbose, "Initial conditions:", part);
    }
    catch (Exception e)
    {
        Log.WriteLine(MessageType.Warning, "Failed to find partition initial conditions, simulation may be unstable.");
        Log.WriteLine(MessageType.Verbose, e.Message);
    }
}

Is there any better way of solving this? It works for most circuits, but I've found one that produces a singular system, not sure why, because it performs fine, event without initial conditions solved 🤔 @dsharlet any ideas?

[dsharlet]

Sorry if this is a dumb suggestion, but can you just provide a better value for the resistances? My inclination would be to take the initial value for the potentiometers (as it will be when simulation begins), and use that to provide the 3 actual values (the wiper variable, and the two resistances corresponding to that wiper value).

If for some reason that is tricky, an easier thing to do might be to simply assume the potentiometers/variable resistors are split in half (so just take the total resistance and divide it by two for the two resistances). For most realistic circuits, I think it probably doesn't matter what the potentiometers are set to for solving for initial conditions (but I could be wrong).

[Federerer]

can you just provide a better value for the resistances?

Sure, but currently solver is not aware of the fact that those are some special unknowns, and I would like to keep it this way, as they might change or we might add other ones and I don't want to have a "dependency" on the solver 😉 My question is, mostly if it's somehow possible to both solve and back substitute just the linear variables in one step without creating another system of equation to perform back substitution, as I have done it. BTW, I've tested a case when I use those calculated conditions as pivot conditions for solving the system later, and it seems to produce a stable result every time (I'm testing random component permutations). Sadly, stable doesn't mean best performing 😒 And it still sometimes produces huge solutions that take ages to compile. Now I need to dig deeper into how SystemOfEquations is implemented to find out WHY.

[dsharlet]

Sure, but currently solver is not aware of the fact that those are some special unknowns, and I would like to keep it this way, as they might change or we might add other ones and I don't want to have a "dependency" on the solver 😉

I guess I'm confused, the resistances associated with potentiometers aren't actually variables solved by the solver during simulation, they're still constants, just symbolic instead of actual numbers. It seems totally reasonable to give them fixed values when solving for the initial conditions?

edit: it seems very similar to what we do with inductors/capacitors (and switches in my PR to hack those).

My question is, mostly if it's somehow possible to both solve and back substitute just the linear variables in one step without creating another system of equation to perform back substitution, as I have done it.

I don't think there is, fusing the two steps into one isn't going to save you anything I don't think.

And it still sometimes produces huge solutions that take ages to compile. Now I need to dig deeper into how SystemOfEquations is implemented to find out WHY.

There may be bugs in this of course, so it's totally worth digging in to if you suspect a problem. That said: this is fundamentally a pretty tricky problem. If you simply write out a symbolic matrix like:

[ a b c x]
[ d e f y ]
[ g h i z ]

and try to row reduce it into an upper triangular matrix, you'll very quickly end up with huge complicated expressions. Back substitution is going to make a bigger mess. A more realistic example for LiveSPICE is something that is mostly constants, but has a few variables scattered around. I think this still quickly turns into a very messy problem. Here's an example with just 3 unknowns (3x3 system) with only 2 symbolic constants (note this is doing both row reduction and back substitution): https://www.wolframalpha.com/input?i=row+reduce+%5B%5B1%2C+a%2C+2%2C+10%5D%2C+%5B3%2C+4%2C+b%2C+11%5D%2C+%5B-2%2C+5%2C+7%2C+12%5D%5D

I think the "size" of those expressions will be O(N^3) with the number of variables in the system (they look a lot like determinants to me, which makes sense as we could get the same result using Cramer's rule), so if you're talking about ~50 variables like a typical LiveSPICE simulation, imagine those expressions being (50/3)^3 = ~5000x bigger! The big wildcard is that maybe the expressions can simplify somewhat. That could be a huge factor in this. I think missed simplifications might be the best hope for this being a tractable problem.

You could try getting the exact system of equations we are trying to solve into wolfram alpha and see what it comes up with, as a point of comparison. I tried doing similar things when I was first trying to get LiveSPICE working. I was trying to check my results against Mathematica, Sage, SymPy, etc., at least at the time, I was struggling to get them to produce results too!

[dsharlet]

I may have accidentally made that example both harder than reality (the matrix isn't sparse) and easier than reality (the variables are all in the strictly upper triangular part of the matrix).

I think in practice, maybe the problem is O(N^2) due to the sparsity of a circuit analysis matrix (not N^3), but that's still going to be a gigantic expression!

[dsharlet]

I found when attempting to update the .net framework version, it is uncovering an interesting bug that might explain some of your struggles with SystemOfEquations:

Unhandled exception: System.Reflection.TargetInvocationException: Exception has been thrown by the target of an invocation.
 ---> System.ArgumentException: Unable to sort because the IComparer.Compare() method returns inconsistent results. Either a value does not compare equal to itself, or one value repeatedly compared to another value yields different results. IComparer: 'System.Comparison`1[System.Int32]'.
   at System.Collections.Generic.ArraySortHelper`1.Sort(Span`1 keys, Comparison`1 comparer)
   at System.MemoryExtensions.Sort[T](Span`1 span, Comparison`1 comparison)
   at System.Linq.EnumerableSorter`2.QuickSort(Int32[] keys, Int32 lo, Int32 hi)
   at System.Linq.EnumerableSorter`1.Sort(TElement[] elements, Int32 count)
   at System.Linq.OrderedEnumerable`1.ToList()
   at ComputerAlgebra.Add.New(IEnumerable`1 Terms) in C:\Users\dshar\GitHub\LiveSPICE\ComputerAlgebra\ComputerAlgebra\Expression\Sum\Add.cs:line 27
   at ComputerAlgebra.SystemOfEquations.Equation.Solve(Expression x) in C:\Users\dshar\GitHub\LiveSPICE\ComputerAlgebra\ComputerAlgebra\Extensions\SystemOfEquations.cs:line 48
   at ComputerAlgebra.SystemOfEquations.Solve(IList`1 x) in C:\Users\dshar\GitHub\LiveSPICE\ComputerAlgebra\ComputerAlgebra\Extensions\SystemOfEquations.cs:line 322
   at Circuit.TransientSolution.Solve(Analysis Analysis, Expression TimeStep, IEnumerable`1 InitialConditions, ILog Log) in C:\Users\dshar\GitHub\LiveSPICE\Circuit\Simulation\TransientSolution.cs:line 0

A bug in our IComparer could explain a lot of weird behavior...

[dsharlet]

I have a fix for this issue in #179. I'd be curious to hear if it resolves some of the issues you were seeing.