Solving System of nonlinear Equations using ModelingToolkit.jl and NonlinearSolve.jl

Specific Domains Numerics
1

Hey there,
I’m fairly new to Julia and this forum (in fact my first post here).
About my problem:
I’m currently trying to solve a system of 3 nonlinear equations that emerged as the solution of the hanging chain problem (following this video).
My studies didn’t involve solving NL systems very much, so I naively defined the system following the MTK documentation. But when I try to solve it I’m only getting a vector of 3 NaNs as solution.
If someone could help and point out my error I would really appreciate it :slight_smile:

In case my versions are relevant:

  • Julia: 1.6.0
  • ModelingToolkit: 5.14.0
  • NonlinearSolve: 0.3.8

Here is my code:

using ModelingToolkit, NonlinearSolve

@variables x y c_1 c_2 λ
@parameters x_1 x_2 y_1 y_2 m g L_0

f(x, y) = c_1/(m*g) * cosh( m*g/c_1*( x + c_2 ) ) - λ/c_1 - y
n       = c_1/(m*g) * ( sinh(m*g/c_1*x_2 + c_2) - sinh(m*g/c_1 * x_1 + c_2) ) - L_0

f_1 = f(x_1, y_1)
f_2 = f(x_2, y_2)

eqs = [
        0 ~ f_1
        0 ~ f_2
        0 ~ n
        ]

ns = NonlinearSystem(eqs, [c_1, c_2, λ], [x_1, x_2, y_1, y_2, m, g, L_0])
guess = [
            c_1 => 1.0
            c_2 => 1.0
            λ   => 1.0
        ]

ps = [
        x_1 => 0.0
        x_2 => 500.0e3
        y_1 => 0.0
        y_2 => 0.0
        m   => 257.0 * 2.7
        g   => 9.81e3
        L_0 => 600.0e3
        ]

prob = NonlinearProblem(ns, guess, ps)
sol = solve(prob, NewtonRaphson())

(Edit before posting this:)
After writing this post it came to my mind, that maybe the System is ill-conditioned with the given parameters. So I set all my current non-zero parameters to 1.0 (L_0 to 1.1 as its the length of the chain) and indeed I got a non-NaN result. Could someone more proficient in the numerics of nonlinear systems point me in the right direction in setting this up please?

Thanks for your time!

1 Like
2

This is a @YingboMa question.

3

Did you get the equations right? I recrafted to using NLSolve on a plain Julia function and also got immediate Infs (see below), but then started looking at the equations, and they don’t make sense dimensionally.

E.g. you have m*g/c_1*( x + c_2 ) as the argument of a cosh and m*g/c_1*x_2 + c_2 as the argument of a sinh. c_2 would have to have dimension of x (length) in the first but be nondimensional in the second. Are the parentheses misplaced or missing?

julia> x1, x2, y1, y2 = 0, 500e03, 0, 0
(0, 500000.0, 0, 0)

julia> m, g, L0 = 257*2.7, 9.81e03, 600e03
(693.9000000000001, 9810.0, 600000.0)

julia> f(c) = [c[1]/(m*g) * cosh(m*g/c[1]*(x1 + c[2])) - c[3]/c[1] - y1; 
               c[1]/(m*g) * cosh(m*g/c[1]*(x2 + c[2])) - c[3]/c[1] - y2; 
               c[1]/(m*g) * (sinh(m*g/c[1]*x2 + c[2]) - sinh(m*g/c[2] * x1 + c[2])) - L0]

f (generic function with 1 method)

julia> using NLsolve

julia> cguess = [1.0; 1.0; 1.0]
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> solution = nlsolve(f, cguess)
ERROR: During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [1, 2, 3]
Stacktrace:
 [1] check_isfinite(x::Vector{Float64})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/utils.jl:24
 [2] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:121
 [3] trust_region (repeats 2 times)
   @ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
 [4] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
 [5] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
 [6] nlsolve(f::Function, initial_x::Vector{Float64})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:46
 [7] top-level scope
   @ REPL[6]:1

julia> f(cguess)
3-element Vector{Float64}:
 Inf
 Inf
 Inf

julia> cguess = [.1; .1; 0.1]
3-element Vector{Float64}:
 0.1
 0.1
 0.1

julia> f(cguess)
3-element Vector{Float64}:
 Inf
 Inf
 Inf
1 Like
4

Thanks, you are right the dimensions of c_2 are wrong.

I think in my function n the term m*g/c_1 (=constant) got absorbed by c_2 and in function f I forgot to do this. I’ll try to fix this and report back.

Edit:
I corrected my f-function, but unfortunately the solution is still a vector of NaNs

5

The issue is that with the choice of units, an order-1 guess for c[1] causes cosh and sinh to overflow. (At first I thought these must be CGS units, but in that case, g = 9.81e02 cm/s^2, right?)

julia> m, g, L0 = 257*2.7, 9.81e03, 600e03;

julia> x1, x2, y1, y2 = 0, 500e03, 0, 0;

julia> c = [1;1;1];

julia> m*g/c[1]*(x1 + c[2])
6.807159000000001e6

julia> cosh(m*g/c[1]*(x1 + c[2]))
Inf

You can prevent the initial overflow by setting the guess for c1 to something large, e.g.

julia> c = [1e06;1;1];

julia> cosh(m*g/c[1]*(x1 + c[2]))
452.14957457816104

But a nonlinear solver is going to have a hard time dealing with a search space where there are large differences of scale on the components of c. Better to nondimensionalize the system or choose units where your numbers are order 1 or order 10. MKS should do.

3 Likes
6

Thanks for pointing out the cosh-function shooting to infinty, guess I lack some intuition with hyperbolas.

The unit system I used was MMGS (realized I made an conversation error in my m term :sweat_smile:).
Switched to MKS for now, but as you said the nonlinear solver doesn’t like the conditioning (is this the right term here?) of the system and throws an LinearAlgebra.SingularException(3) error, even with c_1-guess of 10000. I guess I’ll try to nondimensionalize the system next.

Edit:
I played arround with my c_1-guess and at c_1 = 10 I’m getting a somewhat realistic solution. Thanks again for your help @John_Gibson !

1 Like
7

If you rewrite the argument of cosh and sinh as follows

(m g / c1) (x1 + c2)  = (L m g / c1) (x1/L + c2/L) 

                      =  C1 (x1/L + C2)

with C1 = Lmg/c1 and C2 = c2/L, then C1 and C2 should be order-1 constants. The nonlinear solve over C1, C2 will be better-behaved and will avoid that overflow for order-1 initial guesses.

The right term here is preconditioning: transforming the problem so that it’s better suited to numerical solution methods.

2 Likes
8

Thanks for the clarification.
I’ll try your rewritten arguments out as soon as I find some time for it.