Getting the symbolic solution of system of equations using ModelingToolKit

Hi,

I have the following circuit, which is used to estimate the isolation resistance of high-voltage battery packs to the chassis / low-voltage reference. The switch SW is turned on and off, so there are 2 equations, to solve for RisoP and RisoN, and all the other values are known (the voltage Vadc is sensed by a microncontroller for each state of SW).

I want to get the symbolic equation for bith RisoP and RisoN. I can have the answer using SymPy:

using SymPy

Vpack, Vadc_off, Vadc_on, R1, R2, R3, RisoP, RisoN = symbols("Vpack, Vadc_off, Vadc_on, R1, R2, R3, RisoP, RisoN")

# RisoN || (R2+R3)
A = 1/(1/RisoN + 1/(R2+R3))
# RisoP || R1
B = 1/(1/RisoP + 1/R1)
# R3 / (R2+R3)
Vadc_ratio = R3/(R2+R3)

Vchassis_off = Vpack * (A / (RisoP + A))
Vchassis_on = Vpack * (A / (A+B))

eq1 = Vadc_off ~ Vchassis_off * Vadc_ratio
eq2 = Vadc_on  ~ Vchassis_on * Vadc_ratio

ret = solve((eq1, eq2), (RisoP, RisoN))

RisoP_hat = ret[1][1]
RisoN_hat = ret[1][2]

println(RisoP_hat)
println(RisoN_hat)

----
R1*R3*Vpack*(Vadc_off - Vadc_on)/(Vadc_off*(R2*Vadc_on + R3*Vadc_on - R3*Vpack))
R1*R3*Vpack*(-R2*Vadc_off + R2*Vadc_on - R3*Vadc_off + R3*Vadc_on)/(R1*R3*Vadc_off*Vpack - R1*R3*Vadc_on*Vpack + R2^2*Vadc_off*Vadc_on + 2*R2*R3*Vadc_off*Vadc_on - R2*R3*Vadc_off*Vpack - R2*R3*Vadc_on*Vpack + R3^2*Vadc_off*Vadc_on - R3^2*Vadc_off*Vpack - R3^2*Vadc_on*Vpack + R3^2*Vpack^2)

As I’m trying to use MTK/SciML/Symbolics more extensively, I would like to use it to solve this problem. So, I have 2 questions:

  1. Can SciML’s Symbolics solve that? When I try the following:
using ModelingToolkit
using Groebner

@variables Vpack R1 RisoP RisoN R2 R3 
@variables Vadc_off Vadc_on

# RisoN || (R2+R3)
A = 1/(1/RisoN + 1/(R2+R3))
# RisoP || R1
B = 1/(1/RisoP + 1/R1)
# R3 / (R2+R3)
Vadc_ratio = R3/(R2+R3)

Vchassis_off = Vpack * (A / (RisoP + A))
Vchassis_on = Vpack * (A / (A+B))

eq1 = Vadc_off ~ Vchassis_off * Vadc_ratio
eq2 = Vadc_on  ~ Vchassis_on * Vadc_ratio

symbolic_solve([eq1, eq2], (RisoP, RisoN))

I get

┌ Info: Assuming (R2*RisoN + R2*RisoP + R3*RisoN + R3*RisoP + RisoN*RisoP) != 0
└ @ Symbolics /home/ricardo/.julia/packages/Symbolics/6CYZh/src/solver/preprocess.jl:53
┌ Info: Assuming (R1*R2*RisoN + R1*R2*RisoP + R1*R3*RisoN + R1*R3*RisoP + R1*RisoN*RisoP + R2*RisoN*RisoP + R3*RisoN*RisoP) != 0
└ @ Symbolics /home/ricardo/.julia/packages/Symbolics/6CYZh/src/solver/preprocess.jl:53
"Groebner bases engine is required. Execute `using Groebner` to enable this functionality."
  1. If I model that circuit using @mtkmodel:
using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
using ModelingToolkit: t_nounits as t

@mtkmodel simple_resistor begin
    @parameters begin
        Vpack = 800.0
        R1 = 4.0*499.0e3
        R2 = 4.0*499.0e3
        R3 = 1/3*60.4e3
        RisoP = 50e3
        RisoN = 100e3
    end
    @components begin
        r1    = Resistor(R = R1)
        r2    = Resistor(R = R2)
        r3    = Resistor(R = R3)
        risop = Resistor(R = RisoP)
        rison = Resistor(R = RisoN)
        source = Voltage()
        constant = Constant(k = Vpack)
        ground = Ground()
    end
    @equations begin
        connect(constant.output, source.V)
        connect(source.p, r1.p)
        connect(r1.n, r2.p)
        connect(r2.n, r3.p)

        connect(source.p, risop.p)
        connect(risop.n, rison.p)
        connect(risop.n, r2.p) # connects the resistor network with the chassis

        connect(source.n, r3.n, rison.n, ground.g)
    end
end

@mtkbuild sys_ON = simple_resistor()
@mtkbuild sys_OFF = simple_resistor(R1 = 1e18)

is there a way to get the symbolic solution for a variable, like sys_ON.r3.v and sys_OFF.r3.v (which would be the Vadc when SW is ON and OFF, respecitvely)? With those solutions, I can iterate one more time to solve for RisoP and RisoN.

Did you do this? I assume in your script you did execute that, so did it not load for you?

Yes, I run “using Groebner”, and I can see it’s loaded by:

using Pkg
Pkg.status()
Status `~/.julia/environments/v1.11/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [a93c6f00] DataFrames v1.7.0
⌃ [82cc6244] DataInterpolations v6.5.2
⌃ [0c46a032] DifferentialEquations v7.14.0
  [0b43b601] Groebner v0.8.2
⌃ [7073ff75] IJulia v1.25.0
  [a98d9a8b] Interpolations v0.15.1
⌃ [961ee093] ModelingToolkit v9.49.0
  [16a59e39] ModelingToolkitStandardLibrary v2.17.0
⌅ [8913a72c] NonlinearSolve v3.15.1
⌃ [1dea7af3] OrdinaryDiffEq v6.89.0
⌃ [f0f68f2c] PlotlyJS v0.18.14
⌃ [91a5bcdd] Plots v1.40.8
  [24249f21] SymPy v2.2.0
  [2efcf032] SymbolicIndexingInterface v0.3.34

The full error stack is:

┌ Info: Assuming (R2*RisoN + R2*RisoP + R3*RisoN + R3*RisoP + RisoN*RisoP) != 0
└ @ Symbolics /home/ricardo/.julia/packages/Symbolics/6CYZh/src/solver/preprocess.jl:53
┌ Info: Assuming (R1*R2*RisoN + R1*R2*RisoP + R1*R3*RisoN + R1*R3*RisoP + R1*RisoN*RisoP + R2*RisoN*RisoP + R3*RisoN*RisoP) != 0
└ @ Symbolics /home/ricardo/.julia/packages/Symbolics/6CYZh/src/solver/preprocess.jl:53

"Groebner bases engine is required. Execute `using Groebner` to enable this functionality."

Stacktrace:
 [1] solve_multivar(eqs::Vector{Num}, vars::Tuple{Num, Num}; dropmultiplicity::Bool, warns::Bool)
   @ Symbolics ~/.julia/packages/Symbolics/6CYZh/src/solver/main.jl:328
 [2] symbolic_solve(expr::Vector{Equation}, x::Tuple{Num, Num}; dropmultiplicity::Bool, warns::Bool)
   @ Symbolics ~/.julia/packages/Symbolics/6CYZh/src/solver/main.jl:204
 [3] symbolic_solve(expr::Vector{Equation}, x::Tuple{Num, Num})
   @ Symbolics ~/.julia/packages/Symbolics/6CYZh/src/solver/main.jl:145
 [4] top-level scope
   @ ~/jupyter/isolation_sensing/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W1sZmlsZQ==.jl:1

Open an issue on Symbolics. That is a bug.

I ran into the same problem and opened an issue.

2 Likes

Many thanks for that!

Besides the bug, is it possible to get the analytical solution for a symbol in a @mtkmodel? @ChrisRackauckas

That would be how you do it. And we need to expose more interfaces for better analytical solutions and will do so in the future.

2 Likes