issue with continuation of Neimark-Sacker point

[LEquinoxy]

Hi again!
My problem is related to the ODE in #155: (1) I cannot get the norm form of Neimark-Sacker point (SingualarException will appear); (2) If i passing detect_codim2_bifurcation=2, I cannot continue the Neimark-Sacker point (SingualarException will appear); (3) If i passing detect_codim2_bifurcation=0, the Newton iteration seems extremely slow. What can I do to fix these problems? Please find details in the attachment (the calculation of the Neimark-Sacker point may takes about 10min).
issue with continuation of neimark-sacker point.zip

[rveltz]

are you at JuliaCon by any chance?

[LEquinoxy]

I'm currently in China and unable to attend JuliaCon. 😭 😭 😭 😭

[rveltz]

did you progress on this?

[LEquinoxy]

Sorry for the delayed response! I’ve been working on averaging the original system to obtain a equilibrium-type steady-state for faster continuation, but I'm still curious about why the continuation of the ns point fails. I’ve divided the entire system into two separate parts: the Sending End ODE_LCC_PO and the Receiving End ODE_MMC_PO . The former is significantly stiffer. However, I found that I can calculate the normal form of the Sending End system, but not for the Receiving End system. It seems to get stuck during the initial guess calculation when attempting to continue from the ns point. Please see the details in the attachment. Are there any additional setting that can help the continuation?
BF.zip

[rveltz]

I think your system is very stiff and shooting is not giving you what you think. You should try different solvers at lower tolerance and see if the periodic orbit is correctly captured.

[LEquinoxy]

Thank you for your suggestion! I have one more question, and then perhaps you can close this issue.
From an engineering perspective, certain parameter changes may not lead to significant variations in the periodic orbit. Therefore, some researchers skip the continuation step and directly use the initial periodic orbit to calculate system stability for each parameter. I would like to try using this less precise method of skipping continuation in your BifurcationKit for comparison. Should I just modify the tol in Newton'Par to a large value (like 100)? However, I've noticed that the norm(xtt[:, end],2) still changes for each parameter. Could you provide me with some guidance? Thank you very much!

[rveltz]

This part has piecewise constant part, not sure what newton will give, that's my point.

[rveltz]

I would like to try using this less precise method of skipping continuation in your BifurcationKit for comparison. Should I just modify the

Note sure what you mean. I tend to think that the Shooting continuation of your system is OK (in term of precision) but I question the Floquet exponents. Hence, if this is fast enough for your, I would keep the PALC if you can.

[LEquinoxy]

What I mean is that instead of using PALC(or other methods) to locate the periodic orbit, We
may directly assume that the periodic orbit remains unchanged under parameter variations (i.e., the initial periodic orbit) (in engineering aspect) and then calculate the system's bifurcations. Currently, I'm adjusting tol to a large value so that the first step of the Newton iteration always satisfies the convergence condition for each parameter. The codes related the Norm are shown in Fig 1

Fig 1

The bifurcation digram with tol=100 is shown in Fig 2,

Fig 2

Since the periodic orbit remains unchanged under parameter variations, the corresponding bifurcation diagram should be a straight line, which is clearly different from the results in Fig2. Therefore, simply increasing tol may not be sufficient?

[rveltz]

PALC following the periodic orbits does not seem to be a problem here (so I would not use the constant scheme you are describing). What is more troublesome is the computation of the stability of the periodic orbits. You have to compute the eigenvalues of the derivative of shooting map. However, some components are constant almost everywhere and the precision of the shooting might miss the corresponding eigenvalues. Hence, the detection of bifurcations does not seem very robust for this system and the NS you are seeing is maybe due to lack of precision.

[LEquinoxy]

I see. I will try different solvers and lower tolerances, hoping for some progress😭. You can close this issue for now. Thank you sosososo much!

[LEquinoxy]

Hello again! I recently tried changing the stiff parameter(i.e.,
AA_i and AA_firing), and it seems that the continuation of the NS point can now be done (although setting detect_codim2_bifurcation to 2 still leads to a loop detection of the CH point). That said, plotting the bifurcation curve is already a good start. Now, I’m wondering if it’s possible to use GPU acceleration for the entire computation process (including the continuation of PO and NS). I’ve gone through some of the tutorials on the website, but I didn’t fully understand them. Could you provide more detailed reference materials or specific instructions? Many thanks in advance ☺️

Also, for a periodic system excited by three phase voltage source like this, the frequency of the stable orbit will always be 60(50)Hz, therefore, Is it possible to remove the equation for solving the period? In some cases, if the residual of a Newton iteration step is too large, the computed period in that step might increase, leading to a larger maxiter needed during the shooting process.

[rveltz]

Hi,

Are you still using Shooting nd Matrix-Free? There is very little work in this part for now, I am not sure NS continuation works for matrix-free shooting for now.

I am not sure GPU would lead to acceleration in your case. I would perhaps bet on collocation which is really good now and optimized, you should give it a try.

[LEquinoxy]

Hi, Romain！

Are you still using Shooting nd Matrix-Free?

No, I use jacobian=BK.AutoDiffDense() , since I found that if i pass jacobian=BK.AutoDiffMF() and using PALC, a error:
ERROR: MethodError: no method matching adjoint(::BifurcationKit.var"#800#811"{Vector{…}, @NamedTuple{…}, ShootingProblem{…}}) The function adjoint exists, but no method is defined for this combination of argument types.
will appear at the first newton iteration. If i using Natural(), the newton iteration and continuation will be OK, but the calculated eigenvalue seems wrong.

you should give it a try

I tried the collocation on [email protected], but i will get ERROR: DimensionMismatch: jacobian(f, x) expects that f(x) is an array. Perhaps you meant gradient(f, x)? at the first newton iteration.

Here is my code:

using Revise, Parameters, Plots, BifurcationKit, JLD2,  DifferentialEquations,MATLAB,ParameterizedFunctions,LinearAlgebra

# using Revise, Parameters, Plots, BifurcationKit, JLD2, DifferentialEquations,MATLAB,ParameterizedFunctions,LinearAlgebra

const BK = BifurcationKit
using Arpack
using SparseConnectivityTracer, ADTypes, LinearSolve,TerminalLoggers
#ODE_Setting
RT=1e-8;AT=1e-10
include("ODE_AOC_POrec.jl")
include("plot_sol.jl")


file = jldopen("u0_AOCrec.jld2", "r")
u0 = file["u0_new"]
close(file)
# u0=[zeros(25,1);1;0;zeros(8,1);zeros(38,1)];
# tsidas_alg = VCABM()
tsidas_alg = Vern6()
par_ODE = (AA_i=5e3, AA_firing=5e4, AA_ctrl=5e2, Ihold_pu=2.5e-3, 
           Us_pu_rec=1.0, Kp_dc=1.0989, Ki_dc=1/0.01092, Kp_PLL_rec=10.0, Ki_PLL_rec=1/0.02,
           Us_pu_inv=1.0, Kp_CC=0.5,    Ki_CC=10.0,      Kp_PLL_inv=10.0, Ki_PLL_inv=1/0.02,
           Kp_udc=2.0,Ki_udc=1/0.01,Kp_cQ=0.0,Ki_cQ=1/0.05,
           Lfault_real=0.1,Uref=1.0)

           detector = TracerSparsityDetector()      
           du0 = copy(u0)
           jac_sparsity = ADTypes.jacobian_sparsity(
            (du, u) -> ODE_AOC_POrec(du, u, par_ODE, 0.0), du0, u0, detector)

            condition(u, t, integrator) = u[27]
            affect!(integrator) = nothing
            cb = ContinuousCallback(condition, affect!, nothing,
                save_positions = (true, false))

f = ODEFunction(ODE_AOC_POrec, jac_prototype = float.(jac_sparsity));
prob_ODE = ODEProblem(f, u0, (0,0.5), par_ODE,  reltol=RT, abstol=AT);
sol = @time DifferentialEquations.solve(prob_ODE,progress = true,alg=tsidas_alg)
plot_sol(sol)
u0_new=sol.u[end];
t0_new=sol.t[end];
save("u0_AOCrec.jld2", "u0_new", u0_new)


pinitial=4.0;
Ki_dc_ALL=1/0.01092
par_ODE = @set (par_ODE.Kp_dc=pinitial)
par_ODE = @set (par_ODE.Ki_dc=Ki_dc_ALL[1])

u=hcat(sol.u...)
u0=u[:,end]
 
prob_br = BifurcationProblem(ODE_AOC_POrec, u0, par_ODE, (@optic _.Kp_dc))
prob_ODE = ODEProblem(f, u0, (0.0,0.02), par_ODE,  reltol=RT, abstol=AT);
sol = @time  DifferentialEquations.solve(prob_ODE,progress = true)
##shooting
# probsh, cish = generate_ci_problem(ShootingProblem(M = 1), jacobian=BK.AutoDiffDense(),prob_br, prob_ODE, alg=tsidas_alg,sol, 0.02, reltol=RT, abstol=AT,parallel=true)
# probsh, cish = generate_ci_problem(ShootingProblem(M = 1), jacobian=BK.AutoDiffMF(),prob_br, prob_ODE, alg=tsidas_alg,sol, 0.02, reltol=RT, abstol=AT,parallel=true)    
# @time M1=probsh(cish , par_ODE)##
    # @time Jpo=probsh(Val(:JacobianMatrix), cish,par_ODE)##
    #   using IncompleteLU
    #   Precilu = @time lu(Jpo);
##collocation
probsh, cish =generate_ci_problem( PeriodicOrbitOCollProblem(50, 4; update_section_every_step = 0,meshadapt = true), prob_br, sol, 0.02)

    ls=DefaultLS();       
    # ls = GMRESIterativeSolvers(reltol = 1e-12, N = size(74,1), restart = 20, maxiter =25, Pl = Precilu)
    # ls = GMRESIterativeSolvers(reltol = 1e-7, N = size(74,1), restart = 20, maxiter =25)
	eig = DefaultEig()
	newton_options = NewtonPar(tol =5e-6, max_iterations = 10, linsolver = ls; verbose = true,
	eigsolver=eig,
	linesearch = true,
	)
    opts_br = ContinuationPar(p_min = -50.0, p_max = 20.0, ds = 0.25, dsmax = 0.5, dsmin = 0.25; 
    detect_bifurcation =3,
    n_inversion = 26, 
    max_bisection_steps=15,
    detect_loop=true,
    nev = 3, max_steps = 250, tol_stability = 1e-8, 
    plot_every_step = 1,newton_options)

        argspo = (
            record_from_solution = (x, p; k...) -> begin
		xtt = get_periodic_orbit(p.prob, x, p.p)
		return (max = norm(xtt[:,end],2),
				min = minimum(xtt[1,:]),
				period = getperiod(p.prob, x, p.p))
	end,
	plot_solution = (x, p; k...) -> begin
		xtt = get_periodic_orbit(p.prob, x, p.p)
		plot!(xtt.t, xtt[16,:]; label = "x", k...)
		# plot!(xtt.t, xtt[16,:]; label = "y", k...)
		# plot!(br; subplot = 1, putspecialptlegend = false)
	end)

# alg_cont=DefCont(deflation_operator=DeflationOperator(2, dot, 1.0, [sol.u]),max_branches=1;alg=BK.PALC(),perturb_solution = (x,p,id) -> (x  .+ 0.01 .* rand(length(x))));
# alg_cont=DefCont(deflation_operator=DeflationOperator(2, dot, 1.0, [sol.u]),max_branches=1;alg=BK.Natural(),perturb_solution = (x,p,id) -> (x  .+ 0.01 .* rand(length(x))));
alg_cont=BK.Natural();
# alg_cont=PALC(tangent = Bordered())
# alg_cont=PALC(tangent = Secant())
# alg_cont=MoorePenrose(method = BifurcationKit.iterative)

br2= @time continuation(probsh, cish,alg_cont,
opts_br;
verbosity=3,
plot=true,
bothside=false,
# linear_algo = MatrixBLS(),
callback_newton = (state; linsolver = ls, prob = probsh, p = par_ODE, k...)->begin  if mod(k[:iterationC], 2) == 1 && state.step == 1
                        T=state.x[end];      
                        @info "update period=$T"
                        # p0=k[:z0].p;
                        # @info "parm=$p0"
                        # Jpo = probsh(Val(:JacobianMatrix), state.x,(@set p.Kp_dc = p0))
                        # Precilu = lu(Jpo)
                        # ls.Pl = Precilu
                    end
                    if state.residual > 10
                    R=state.residual
                    @warn "Reject Newton step!!(the residual is too large: residual=$R)"
                    return false
                    end

                    if state.x[end]>0.021
                    T=state.x[end];
                    @warn "Reject Newton step!!(the period is too large: period=$T)"
                    return false  
                    end

                    return true
                    end,

finalise_solution = (z, tau, step, contResult; k...) -> begin
return !any(sp -> sp.ind_ev >= 2, contResult.specialpoint)
end,
normC = norminf,
argspo...)

TEST1_0116.zip