ODE solving is fast on DifferentialEquation@7.9.0 but significantly slower on later versions

Specific Domains Numerics
,
1

When solving the same ODE (a 100-dimensional strongly slow-fast system, where the non-smooth component is approximated using hyperbolic tangent function ) problem using DifferentialEquations.jl , I noticed that ODE solving is fast on DifferentialEquation@7.9.0 but significantly slower on later versions

DifferentialEquation@7.9.0:

BenchmarkTools.Trial: 10 samples with 1 evaluation.
 Range (min … max):  1.883 s …    2.348 s  ┊ GC (min … max):  7.62% … 22.09%
 Time  (median):     2.093 s               ┊ GC (median):     9.50%
 Time  (mean ± σ):   2.118 s ± 156.198 ms  ┊ GC (mean ± σ):  11.74% ±  8.24%

  █    █            ██  █      █          █     █    █     █
  █▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁██▁▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁█▁▁▁▁▁█ ▁
  1.88 s         Histogram: frequency by time         2.35 s <

 Memory estimate: 595.07 MiB, allocs estimate: 2165905.

DifferentialEquation@7.14.0:

BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max):  16.519 s …   16.991 s  ┊ GC (min … max): 2.12% … 4.28%
 Time  (median):     16.755 s               ┊ GC (median):    3.21%
 Time  (mean ± σ):   16.755 s ± 333.446 ms  ┊ GC (mean ± σ):  3.21% ± 1.53%

  █                                                         █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  16.5 s          Histogram: frequency by time           17 s <

 Memory estimate: 482.70 MiB, allocs estimate: 3585645.

The julia version info:

Julia Version 1.10.4
Commit 48d4fd4843 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 24 × Intel(R) Core(TM) i7-14650HX
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 24 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =

The code:

using JLD2,DifferentialEquations
#ODE_Setting
RT=1e-8;AT=1e-10

include("ODE_SOC_PO.jl")

using BenchmarkTools

 u0=[zeros(25,1);1;0;zeros(8,1);zeros(38,1)];

tsidas_alg = AutoVern6(Rodas5P(),nonstifftol=15/10)
# AutoVern7(ODEInterfaceDiffEq.Rodas5(),nonstifftol = 15/ 10)# Decreasing nonstifftol makes switching to the stiff algorithm more likely
par_ODE = (AA_i=5e4, AA_firing=5e5, AA_ctrl=5e2, Ihold_pu=2.5e-4, 
           Us_pu_rec=0.82, 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)
prob_ODE = ODEProblem(ODE_SOC_PO, u0, (3,4), par_ODE, reltol=RT, abstol=AT)
a=@time DifferentialEquations.solve(prob_ODE,progress = true,alg=tsidas_alg)

The ODE:

function ODE_SOC_PO(dydt, X, p, t = 0)
	AA_i, AA_firing, AA_ctrl, Ihold_pu, 
    Us_pu_rec, Kp_dc, Ki_dc, Kp_PLL_rec, Ki_PLL_rec,
    Us_pu_inv, Kp_CC, Ki_CC, Kp_PLL_inv, Ki_PLL_inv,
    Kp_udc,Ki_udc,Kp_cQ,Ki_cQ,
    Lfault_real = p

Uref=1.0;

idrefmax=1.5;
idrefmin=-1.5;

iqrefmax=1.5;
iqrefmin=-1.5;
del=1e-5;

cof=1e-5;
# cof=1e-3;

    y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,
    y11,y12,y13,y14,y15,y16,y17,y18,y19,y20,
    y21,y22,y23,y24,y25,y26,y27,y28,y29,y30,
    y31,y32,y33,y34,y35,y36,y37,y38,y39,y40,
    y41,y42,y43,y44,y45,y46,y47,y48,y49,y50,
    y51,y52,y53,y54,y55,y56,y57,y58,y59,y60,
    y61,y62,y63,y64,y65,y66,y67,y68,y69,y70,
    y71,y72,y73=X
    

    COS=y26;
    SIN=y27;
    COS_N=y26;
    SIN_N=-y27;
    COS_2=1-2*y27^2;
    SIN_2=2*y26*y27;
    COS_2N=1-2*y27^2;
    SIN_2N=-2*y26*y27;

    w0=100*pi;
    Sbase_rec=603.73;
    Ubase1_L2L_rec=345;
    Ubase2_L2L_rec=213.457;
    Ubase1_P2P_rec=Ubase1_L2L_rec*sqrt(2/3);
    Ubase2_P2P_rec=Ubase2_L2L_rec*sqrt(2/3);
    Ibase2_P2P_rec=Sbase_rec/Ubase2_L2L_rec*sqrt(2/3);
    Ibase1_P2P_rec=Sbase_rec/Ubase1_L2L_rec*sqrt(2/3);
    Kt_rec=Ubase1_L2L_rec/Ubase2_L2L_rec;Tra_rec=Ubase2_L2L_rec^2/Ubase1_L2L_rec^2;
    #####
    Sbase_inv=1440;Ubase1_L2L_inv=525;Ubase2_L2L_inv=244;
    Ubase2_P2P_inv=Ubase2_L2L_inv*sqrt(2/3);
    Ibase2_P2P_inv=Sbase_inv/Ubase2_L2L_inv*sqrt(2/3);
    Kt_inv=Ubase1_L2L_inv/Ubase2_L2L_inv;Tra_inv=Ubase2_L2L_inv^2/Ubase1_L2L_inv^2;
    #####

    Lt_rec=0.18*Ubase2_L2L_rec^2/Sbase_rec/w0/(Ubase2_L2L_rec^2/Sbase_rec);
    Ihold=Ihold_pu;
    Crc=0.0175e-6*(Ubase2_L2L_rec^2/Sbase_rec);
    Rrc=2880.0/(Ubase2_L2L_rec^2/Sbase_rec);
    # Lg_rec=0.07*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    # Rg_rec=2160.633*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Ldc_rec=0.3/(Ubase2_L2L_rec^2/Sbase_rec);

    Lt_inv=0.18/w0;
    Larm=0.04/(Ubase2_L2L_inv^2/Sbase_inv);
    Ceq=(18000e-6)/200*(Ubase2_L2L_inv^2/Sbase_inv);
    Leq=Lt_inv+0.5*Larm;
    Kd=0.5*100*pi*Larm+0.18;
    Lg_inv=0.037*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    Rg_inv=45.11*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    Rg_inv_2=2.912*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    Rdc_inv=3/(Ubase2_L2L_inv^2/Sbase_inv);
    Ldc_inv=0.15/(Ubase2_L2L_inv^2/Sbase_inv);
    Cdc_inv=1e-6*(Ubase2_L2L_inv^2/Sbase_inv);

    Rline=2.5/(Ubase2_L2L_rec^2/Sbase_rec);
    Cline=26e-6*(Ubase2_L2L_rec^2/Sbase_rec);

    Usa_rec=Us_pu_rec*525*sqrt(2/3)*COS/(Ubase2_L2L_rec*sqrt(2/3));
    Usb_rec=Us_pu_rec*525*sqrt(2/3)*(COS*cos(pi*2/3)+SIN*sin(pi*2/3))/(Ubase2_L2L_rec*sqrt(2/3));
    Usc_rec=-Usa_rec-Usb_rec;
    Usa_inv=Us_pu_inv*525*sqrt(2/3)*COS/Ubase2_P2P_inv;
    Usb_inv=Us_pu_inv*525*sqrt(2/3)*(COS*cos(-2*pi/3)-SIN*sin(-2*pi/3))/Ubase2_P2P_inv;
    Usc_inv=-Usa_inv-Usb_inv;

    ###REC
    idc_rec_end=y25;
    Edc_rec=y70;
    Upcca_rec=y13;
    Upccb_rec=y14;
    Upccc_rec=y34;
    ia_Y=y19;ib_Y=y20;ic_Y=-y19-y20;
    ia_D=y22;ib_D=y23;ic_D=-y22-y23;
    ##INV
    idc_inv_end=Ibase2_P2P_inv*y69/Ibase2_P2P_rec
    Edc_inv=Ubase2_P2P_rec*(Edc_rec+Rline*idc_inv_end)/(Ubase2_P2P_inv);

    # Upcca_inv=Usa_inv/Kt_inv-Rg_inv*((y36+y73)-y47);
    # Upccb_inv=Usb_inv/Kt_inv-Rg_inv*(y37-y48);
    # Upccc_inv=Usc_inv/Kt_inv-Rg_inv*((-y36-y37)-y72);

    # Upcca_inv=Usa_inv/Kt_inv-Rg_inv*(y36-y47);
    # Upccb_inv=Usb_inv/Kt_inv-Rg_inv*(y37-y48);
    # Upccc_inv=-Upcca_inv-Upccb_inv;

    Upcca_inv=Usa_inv/Kt_inv-(Rg_inv_2*(y36+y73)+Rg_inv*((y36+y73)-y47));
    Upccb_inv=Usb_inv/Kt_inv-(Rg_inv_2*(y37)+Rg_inv*(y37-y48));
    Upccc_inv=Usc_inv/Kt_inv-(Rg_inv_2*(-y36-y37)+Rg_inv*((-y36-y37)-y72));

    udc_inv=y68-Rdc_inv*(y38+y39+y40);
    T=y58;

    Rg1_rec=3.737*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Rg2_rec=2160.633*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Lg_rec=0.151*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    dydt[1]=1/(Lg_rec+Rg1_rec/Rg2_rec*Lg_rec)*(-Rg1_rec*y1+Usa_rec/Kt_rec-Upcca_rec);
    dydt[2]=1/(Lg_rec+Rg1_rec/Rg2_rec*Lg_rec)*(-Rg1_rec*y2+Usb_rec/Kt_rec-Upccb_rec);
    dydt[28]=1/(Lg_rec+Rg1_rec/Rg2_rec*Lg_rec)*(-Rg1_rec*y28+Usc_rec/Kt_rec-Upccc_rec);

    Rf1=83.32*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Cf1=6.685e-6/Tra_rec*(Ubase2_L2L_rec^2/Sbase_rec);
    Lf1=0.0136*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    dydt[3]=1/Lf1*(Upcca_rec-y5);
    dydt[4]=1/Lf1*(Upccb_rec-y6);
    dydt[29]=1/Lf1*(Upccc_rec-y30);
    dydt[5]=1/Cf1*(y3+1/Rf1*(Upcca_rec-y5));
    dydt[6]=1/Cf1*(y4+1/Rf1*(Upccb_rec-y6));
    dydt[30]=1/Cf1*(y29+1/Rf1*(Upccc_rec-y30));

    Rf21=261.87*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Rf22=29.76*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Lf2=0.1364*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    Cf21=6.685e-6/Tra_rec*(Ubase2_L2L_rec^2/Sbase_rec);
    Cf22=74.28e-6/Tra_rec*(Ubase2_L2L_rec^2/Sbase_rec);
    dydt[7]=1/Cf21*(y11+1/Rf21*(Upcca_rec-y7));
    dydt[8]=1/Cf21*(y12+1/Rf21*(Upccb_rec-y8));
    dydt[31]=1/Cf21*(y33+1/Rf21*(Upccc_rec-y31));
    dydt[9]=1/Cf22*y11;
    dydt[10]=1/Cf22*y12;
    dydt[32]=1/Cf22*y33;
    dydt[11]=1/Lf2*(Upcca_rec-Rf22*y11-(y7+y9));
    dydt[12]=1/Lf2*(Upccb_rec-Rf22*y12-(y8+y10));
    dydt[33]=1/Lf2*(Upccc_rec-Rf22*y33-(y31+y32));

    num_C=1;
    Cf3=3.342e-6/Tra_rec*(Ubase2_L2L_rec^2/Sbase_rec)*num_C;
    ima=(Rg2_rec*y1 +(Usa_rec/Kt_rec-Upcca_rec))/(Rg1_rec+Rg2_rec)-(y3 +1/Rf1*(Upcca_rec-y5 ))-(y11+1/Rf21*(Upcca_rec-y7 ))-(ia_Y+( (3^(1/2)*ia_D)/3 - (3^(1/2)*ib_D)/3))-y35;
    imb=(Rg2_rec*y2 +(Usb_rec/Kt_rec-Upccb_rec))/(Rg1_rec+Rg2_rec)-(y4 +1/Rf1*(Upccb_rec-y6 ))-(y12+1/Rf21*(Upccb_rec-y8 ))-(ib_Y+( (3^(1/2)*ib_D)/3 - (3^(1/2)*ic_D)/3));
    imc=(Rg2_rec*y28+(Usc_rec/Kt_rec-Upccc_rec))/(Rg1_rec+Rg2_rec)-(y29+1/Rf1*(Upccc_rec-y30))-(y33+1/Rf21*(Upccc_rec-y31))-(ic_Y+(-(3^(1/2)*ia_D)/3 + (3^(1/2)*ic_D)/3));
    dydt[13]=1/Cf3*(ima);
    dydt[14]=1/Cf3*(imb);
    dydt[34]=1/Cf3*(imc);
    ######
    Lfault=Lfault_real*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    # Lfault_real*Tra_rec/(Ubase2_L2L_rec^2/Sbase_rec);
    dydt[35]=-10*y35;
    # dydt[35]=1/Lfault*Upcca_rec;
##VDCOL
    Gvdcol=0.0023;
    Tvdcol=0.03;
    dydt[71]=1/Tvdcol*(Gvdcol*Ubase2_P2P_rec*y70-y71);
    svdcol_1=0.5*(1+tanh(-AA_ctrl*(y71)));svdcol_2=0.5*(1+tanh(AA_ctrl*(y71-0.95)));
    Iref=svdcol_1*0.1+svdcol_2*1+(1-(svdcol_1+svdcol_2)).*(18/19*y71+0.1);
    # Iref=1;

    # Tm=0.004;
    # Gm=0.32;
    # dydt[15]=-1/Tm*y15+Gm/Tm*(Ibase2_P2P_rec*y25);
    # dydt[16]=Ki_dc*(Iref-y15);
    # Alpha_ori=pi-(Kp_dc*(Iref-y15)+y16);

    # ind_alpha_Upper=((1+tanh(AA_ctrl*(Alpha_ori-3.054)))/2);
    # ind_alpha_Lower=((1+tanh(AA_ctrl*(0.1-Alpha_ori)))/2);
    # Alpha=3.054*ind_alpha_Upper+0.1*ind_alpha_Lower+(1-(ind_alpha_Upper+ind_alpha_Lower))*Alpha_ori;

    # PI_rec_max=pi-0.1;
    # PI_rec_min=pi-3.054;
    # Tm=0.001;
    # Gm=0.32;
    # AA_ctrl_2=1e4;
    # dydt[15]=-1/Tm*y15+Gm/Tm*(Ibase2_P2P_rec*y25);
    # Alpha_ori=(Kp_dc*(Iref-y15)+y16);
    
    # ind_alp_Up=((1+tanh(AA_ctrl*(Alpha_ori-PI_rec_max)))/2);
    # ind_alp_Lo=((1+tanh(AA_ctrl*(PI_rec_min-Alpha_ori)))/2);

    # ind_alp_Up_inn=((1+tanh(AA_ctrl_2*(y16-(3.054-del))))/2);
    # ind_alp_Lo_inn=((1+tanh(AA_ctrl_2*((0.1+del)-y16)))/2);

    # Lim_alpha=-10*(y16-3)*ind_alp_Up_inn-10*(y16-0.1)*ind_alp_Lo_inn+(1-(ind_alp_Up_inn+ind_alp_Lo_inn))*Ki_dc*(Iref-y15)

    # dydt[16]=(1-(ind_alp_Up+ind_alp_Lo))*Ki_dc*(Iref-y15)+(ind_alp_Up+ind_alp_Lo)*(0);
    # # dydt[16]=Lim_alpha;
    
    # # (1-(ind_alp_Up+ind_alp_Lo))*Ki_dc*(Iref-y15)+(ind_alp_Up+ind_alp_Lo)*(0);
    # Alpha=3.054*ind_alp_Up+0.1*ind_alp_Lo+(1-(ind_alp_Up+ind_alp_Lo))*Alpha_ori;


    PI_rec_max=pi-0.1;
    PI_rec_min=pi-(pi/2-0.01);
    # PI_rec_min=pi-(3.054);
    Tm=0.004;
    Gm=0.32;
    dydt[15]=-1/Tm*y15+Gm/Tm*(Ibase2_P2P_rec*y25);
    Alpha_ori=Kp_dc*(Iref-y15)+y16;
     
    ind_alp_Up=((1+tanh(AA_ctrl*(Alpha_ori-PI_rec_max)))/2);
    ind_alp_Lo=((1+tanh(AA_ctrl*(PI_rec_min-Alpha_ori)))/2);
    Ki_dc_rea=(1-(ind_alp_Up+ind_alp_Lo))*Ki_dc+(ind_alp_Up+ind_alp_Lo)*cof;
    # Ki_dc_rea=Ki_dc;
    dydt[16]=Ki_dc_rea*(Iref-y15);
    # dydt[16]=Ki_dc*(Iref-y15);
    Alpha=pi-(PI_rec_max*ind_alp_Up+PI_rec_min*ind_alp_Lo+(1-(ind_alp_Up+ind_alp_Lo))*Alpha_ori);

    dydt[17]=Ubase1_P2P_rec*Ki_PLL_rec*2/3*(-(sin(y18)*COS+cos(y18)*SIN)*y13-(sin(y18-2/3*pi)*COS+cos(y18-2/3*pi)*SIN)*y14-(sin(y18+2/3*pi)*COS+cos(y18+2/3*pi)*SIN)*(y34));
    dydt[18]=y17+Ubase1_P2P_rec*Kp_PLL_rec*2/3*(-(sin(y18)*COS+cos(y18)*SIN)*y13-(sin(y18-2/3*pi)*COS+cos(y18-2/3*pi)*SIN)*y14-(sin(y18+2/3*pi)*COS+cos(y18+2/3*pi)*SIN)*(y34));
    Pshift1_Y=y18+pi/3;Pshift2_Y=y18;Pshift3_Y=y18-pi/3;Pshift4_Y=y18-2*pi/3;Pshift5_Y=y18-pi;Pshift6_Y=y18-4*pi/3;
    Pshift1_D=y18+pi/3-pi/6;Pshift2_D=y18-pi/6;Pshift3_D=y18-pi/3-pi/6;Pshift4_D=y18-2*pi/3-pi/6;Pshift5_D=y18-pi-pi/6;Pshift6_D=y18-4*pi/3-pi/6;

    Alpha2=Alpha+pi/6; 

    #########################
    C11_Y=COS*cos(Pshift1_Y-Alpha/2)-SIN*sin(Pshift1_Y-Alpha/2);C12_Y=COS*cos(Pshift1_Y-Alpha2/2)-SIN*sin(Pshift1_Y-Alpha2/2);
    Su1_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C11_Y)));Su1_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C12_Y)));
    Su1_Y=Su1_Y_part1-Su1_Y_part2;
    #########################
    C21_Y=COS*cos(Pshift2_Y-Alpha/2)-SIN*sin(Pshift2_Y-Alpha/2);C22_Y=COS*cos(Pshift2_Y-Alpha2/2)-SIN*sin(Pshift2_Y-Alpha2/2);
    Su2_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C21_Y)));Su2_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C22_Y)));
    Su2_Y=Su2_Y_part1-Su2_Y_part2;
    #########################
    C31_Y=COS*cos(Pshift3_Y-Alpha/2)-SIN*sin(Pshift3_Y-Alpha/2);C32_Y=COS*cos(Pshift3_Y-Alpha2/2)-SIN*sin(Pshift3_Y-Alpha2/2);
    Su3_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C31_Y)));Su3_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C32_Y)));
    Su3_Y=Su3_Y_part1-Su3_Y_part2;
    #########################
    C41_Y=COS*cos(Pshift4_Y-Alpha/2)-SIN*sin(Pshift4_Y-Alpha/2);C42_Y=COS*cos(Pshift4_Y-Alpha2/2)-SIN*sin(Pshift4_Y-Alpha2/2);
    Su4_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C41_Y)));Su4_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C42_Y)));
    Su4_Y=Su4_Y_part1-Su4_Y_part2;
    #########################
    C51_Y=COS*cos(Pshift5_Y-Alpha/2)-SIN*sin(Pshift5_Y-Alpha/2);C52_Y=COS*cos(Pshift5_Y-Alpha2/2)-SIN*sin(Pshift5_Y-Alpha2/2);
    Su5_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C51_Y)));Su5_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C52_Y)));
    Su5_Y=Su5_Y_part1-Su5_Y_part2;
    #########################
    C61_Y=COS*cos(Pshift6_Y-Alpha/2)-SIN*sin(Pshift6_Y-Alpha/2);C62_Y=COS*cos(Pshift6_Y-Alpha2/2)-SIN*sin(Pshift6_Y-Alpha2/2);
    Su6_Y_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C61_Y)));Su6_Y_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C62_Y)));
    Su6_Y=Su6_Y_part1-Su6_Y_part2;
    #########################
    Si1_Y = 0.5*(1+tanh(AA_i*(ia_Y-Ihold)));Si3_Y = 0.5*(1+tanh(AA_i*(ib_Y-Ihold)));Si5_Y = 0.5*(1+tanh(AA_i*(ic_Y-Ihold)));
    Si4_Y = 0.5*(1+tanh(AA_i*(-ia_Y-Ihold)));Si6_Y = 0.5*(1+tanh(AA_i*(-ib_Y-Ihold)));Si2_Y = 0.5*(1+tanh(AA_i*(-ic_Y-Ihold)));
    #########################
    S1_Y=Su1_Y+Si1_Y-Su1_Y*Si1_Y;S2_Y=Su2_Y+Si2_Y-Su2_Y*Si2_Y;S3_Y=Su3_Y+Si3_Y-Su3_Y*Si3_Y;
    S4_Y=Su4_Y+Si4_Y-Su4_Y*Si4_Y;S5_Y=Su5_Y+Si5_Y-Su5_Y*Si5_Y;S6_Y=Su6_Y+Si6_Y-Su6_Y*Si6_Y;
    Sa_tri_Y=S1_Y-S4_Y;Sb_tri_Y=S3_Y-S6_Y;Sc_tri_Y=S5_Y-S2_Y;
    # Sa_mul_Y=Sa_tri_Y^2;Sb_mul_Y=Sb_tri_Y^2;Sc_mul_Y=Sc_tri_Y^2;
    Sa_mul_Y=S1_Y+S4_Y-S1_Y*S4_Y;
    Sb_mul_Y=S3_Y+S6_Y-S3_Y*S6_Y;
    Sc_mul_Y=S5_Y+S2_Y-S5_Y*S2_Y;
    #########################
    udc_Y=(Rrc*(0.5*(Sa_tri_Y*ia_Y+Sb_tri_Y*ib_Y+Sc_tri_Y*ic_Y)-y25)+y21);

    ua_Y=Upcca_rec;
    ub_Y=Upccb_rec;
    uc_Y=Upccc_rec;
 
  
    dydt[19]=-(udc_Y*(Sa_tri_Y/2 - (Sa_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/6) - (Sa_mul_Y*ua_Y)/2 + (Sa_mul_Y*Sb_mul_Y*ub_Y)/2 + (Sa_mul_Y*Sc_mul_Y*uc_Y)/2)/Lt_rec;
    dydt[20]=-(udc_Y*(Sb_tri_Y/2 - (Sb_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/6) - (Sb_mul_Y*ub_Y)/2 + (Sa_mul_Y*Sb_mul_Y*ua_Y)/2 + (Sb_mul_Y*Sc_mul_Y*uc_Y)/2)/Lt_rec;
  
    # dydt[19]=-(udc_Y*(Sa_tri_Y/2 - (Sa_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/(2*(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))) + (Sa_mul_Y*Sb_mul_Y*ub_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) + (Sa_mul_Y*Sc_mul_Y*uc_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) - (Sa_mul_Y*ua_Y*(Sb_mul_Y + Sc_mul_Y))/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))/Lt_rec;
    # dydt[20]=-(udc_Y*(Sb_tri_Y/2 - (Sb_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/(2*(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))) + (Sa_mul_Y*Sb_mul_Y*ua_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) + (Sb_mul_Y*Sc_mul_Y*uc_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) - (Sb_mul_Y*ub_Y*(Sa_mul_Y + Sc_mul_Y))/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))/Lt_rec;

    # dydt[19]=-(udc_Y*(Sa_tri_Y/2 - (Sa_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/(2*(3))) + (Sa_mul_Y*Sb_mul_Y*ub_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) + (Sa_mul_Y*Sc_mul_Y*uc_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) - (Sa_mul_Y*ua_Y*(Sb_mul_Y + Sc_mul_Y))/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))/Lt_rec;
    # dydt[20]=-(udc_Y*(Sb_tri_Y/2 - (Sb_mul_Y*(Sa_tri_Y + Sb_tri_Y + Sc_tri_Y))/(2*(3))) + (Sa_mul_Y*Sb_mul_Y*ua_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) + (Sb_mul_Y*Sc_mul_Y*uc_Y)/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y) - (Sb_mul_Y*ub_Y*(Sa_mul_Y + Sc_mul_Y))/(Sa_mul_Y + Sb_mul_Y + Sc_mul_Y))/Lt_rec;
    dydt[21]=1/Crc*(0.5*(Sa_tri_Y*ia_Y+Sb_tri_Y*ib_Y+Sc_tri_Y*ic_Y)-y25);

    C11_D=COS*cos(Pshift1_D-Alpha/2)-SIN*sin(Pshift1_D-Alpha/2);C12_D=COS*cos(Pshift1_D-Alpha2/2)-SIN*sin(Pshift1_D-Alpha2/2);
    Su1_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C11_D)));Su1_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C12_D)));
    Su1_D=Su1_D_part1-Su1_D_part2;
    #########################
    C21_D=COS*cos(Pshift2_D-Alpha/2)-SIN*sin(Pshift2_D-Alpha/2);C22_D=COS*cos(Pshift2_D-Alpha2/2)-SIN*sin(Pshift2_D-Alpha2/2);
    Su2_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C21_D)));Su2_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C22_D)));
    Su2_D=Su2_D_part1-Su2_D_part2;
    #########################
    C31_D=COS*cos(Pshift3_D-Alpha/2)-SIN*sin(Pshift3_D-Alpha/2);C32_D=COS*cos(Pshift3_D-Alpha2/2)-SIN*sin(Pshift3_D-Alpha2/2);
    Su3_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C31_D)));Su3_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C32_D)));
    Su3_D=Su3_D_part1-Su3_D_part2;
    #########################
    C41_D=COS*cos(Pshift4_D-Alpha/2)-SIN*sin(Pshift4_D-Alpha/2);C42_D=COS*cos(Pshift4_D-Alpha2/2)-SIN*sin(Pshift4_D-Alpha2/2);
    Su4_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C41_D)));Su4_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C42_D)));
    Su4_D=Su4_D_part1-Su4_D_part2;
    #########################
    C51_D=COS*cos(Pshift5_D-Alpha/2)-SIN*sin(Pshift5_D-Alpha/2);C52_D=COS*cos(Pshift5_D-Alpha2/2)-SIN*sin(Pshift5_D-Alpha2/2);
    Su5_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C51_D)));Su5_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C52_D)));
    Su5_D=Su5_D_part1-Su5_D_part2;
    #########################
    C61_D=COS*cos(Pshift6_D-Alpha/2)-SIN*sin(Pshift6_D-Alpha/2);C62_D=COS*cos(Pshift6_D-Alpha2/2)-SIN*sin(Pshift6_D-Alpha2/2);
    Su6_D_part1=0.5*(1+tanh(AA_firing*(cos(Alpha/2)-C61_D)));Su6_D_part2=0.5*(1+tanh(AA_firing*(cos(Alpha2/2)-C62_D)));
    Su6_D=Su6_D_part1-Su6_D_part2;
    #########################
    Si1_D=0.5*(1+tanh(AA_i*(ia_D-Ihold)));Si3_D=0.5*(1+tanh(AA_i*(ib_D-Ihold)));Si5_D=0.5*(1+tanh(AA_i*(ic_D-Ihold)));
    Si4_D=0.5*(1+tanh(AA_i*(-ia_D-Ihold)));Si6_D=0.5*(1+tanh(AA_i*(-ib_D-Ihold)));Si2_D=0.5*(1+tanh(AA_i*(-ic_D-Ihold)));
    #########################
    S1_D=Su1_D+Si1_D-Su1_D*Si1_D;S2_D=Su2_D+Si2_D-Su2_D*Si2_D;S3_D=Su3_D+Si3_D-Su3_D*Si3_D;
    S4_D=Su4_D+Si4_D-Su4_D*Si4_D;S5_D=Su5_D+Si5_D-Su5_D*Si5_D;S6_D=Su6_D+Si6_D-Su6_D*Si6_D;
    Sa_tri_D=S1_D-S4_D;Sb_tri_D=S3_D-S6_D;Sc_tri_D=S5_D-S2_D;
    # Sa_mul_D=Sa_tri_D^2;Sb_mul_D=Sb_tri_D^2;Sc_mul_D=Sc_tri_D^2;
    Sa_mul_D=S1_D+S4_D-S1_D*S4_D;
    Sb_mul_D=S3_D+S6_D-S3_D*S6_D;
    Sc_mul_D=S5_D+S2_D-S5_D*S2_D;
    #########################
    udc_D=(Rrc*(0.5*(Sa_tri_D*ia_D+Sb_tri_D*ib_D+Sc_tri_D*ic_D)-y25)+y24);

    ua_D=Upcca_rec/3 + Upccb_rec/3 + Upccc_rec/3 + (3^(1/2)*Upcca_rec)/3 - (3^(1/2)*Upccc_rec)/3;
    ub_D=Upcca_rec/3 + Upccb_rec/3 + Upccc_rec/3 - (3^(1/2)*Upcca_rec)/3 + (3^(1/2)*Upccb_rec)/3;
    uc_D=Upcca_rec/3 + Upccb_rec/3 + Upccc_rec/3 - (3^(1/2)*Upccb_rec)/3 + (3^(1/2)*Upccc_rec)/3;

    dydt[22]=-(udc_D*(Sa_tri_D/2 - (Sa_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/6) - (Sa_mul_D*ua_D)/2 + (Sa_mul_D*Sb_mul_D*ub_D)/2 + (Sa_mul_D*Sc_mul_D*uc_D)/2)/Lt_rec;
    dydt[23]=-(udc_D*(Sb_tri_D/2 - (Sb_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/6) - (Sb_mul_D*ub_D)/2 + (Sa_mul_D*Sb_mul_D*ua_D)/2 + (Sb_mul_D*Sc_mul_D*uc_D)/2)/Lt_rec;

    # dydt[22]=-(udc_D*(Sa_tri_D/2 - (Sa_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/(2*(Sa_mul_D + Sb_mul_D + Sc_mul_D))) + (Sa_mul_D*Sb_mul_D*ub_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) + (Sa_mul_D*Sc_mul_D*uc_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) - (Sa_mul_D*ua_D*(Sb_mul_D + Sc_mul_D))/(Sa_mul_D + Sb_mul_D + Sc_mul_D))/Lt_rec;
    # dydt[23]=-(udc_D*(Sb_tri_D/2 - (Sb_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/(2*(Sa_mul_D + Sb_mul_D + Sc_mul_D))) + (Sa_mul_D*Sb_mul_D*ua_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) + (Sb_mul_D*Sc_mul_D*uc_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) - (Sb_mul_D*ub_D*(Sa_mul_D + Sc_mul_D))/(Sa_mul_D + Sb_mul_D + Sc_mul_D))/Lt_rec;

    # dydt[22]=-(udc_D*(Sa_tri_D/2 - (Sa_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/(2*(3))) + (Sa_mul_D*Sb_mul_D*ub_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) + (Sa_mul_D*Sc_mul_D*uc_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) - (Sa_mul_D*ua_D*(Sb_mul_D + Sc_mul_D))/(Sa_mul_D + Sb_mul_D + Sc_mul_D))/Lt_rec;
    # dydt[23]=-(udc_D*(Sb_tri_D/2 - (Sb_mul_D*(Sa_tri_D + Sb_tri_D + Sc_tri_D))/(2*(3))) + (Sa_mul_D*Sb_mul_D*ua_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) + (Sb_mul_D*Sc_mul_D*uc_D)/(Sa_mul_D + Sb_mul_D + Sc_mul_D) - (Sb_mul_D*ub_D*(Sa_mul_D + Sc_mul_D))/(Sa_mul_D + Sb_mul_D + Sc_mul_D))/Lt_rec;
    dydt[24]=1/Crc*(0.5*(Sa_tri_D*ia_D+Sb_tri_D*ib_D+Sc_tri_D*ic_D)-y25);

    dydt[25]=1/Ldc_rec*((udc_Y+udc_D)-(Edc_rec+Rline*idc_rec_end));
#### COS及SIN
    dydt[26]=-w0*y27+y26*(1-y26^2-y27^2);
    dydt[27]= w0*y26+y27*(1-y26^2-y27^2);

    ## Inverter side ODE
    #PCC voltage/current in positive/negative dq frame (inner variables)
    #idq_P_pu
    id_P_pu= 2/3*((cos(T)*COS-sin(T)*SIN)*y36+(cos(T-2/3*pi)*COS-sin(T-2/3*pi)*SIN)*y37+(cos(T+2/3*pi)*COS-sin(T+2/3*pi)*SIN)*(-y36-y37));
    iq_P_pu=-2/3*((sin(T)*COS+cos(T)*SIN)*y36+(sin(T-2/3*pi)*COS+cos(T-2/3*pi)*SIN)*y37+(sin(T+2/3*pi)*COS+cos(T+2/3*pi)*SIN)*(-y36-y37));
    #udq_P_pu
    ud_P_pu= 2/3*((cos(T)*COS-sin(T)*SIN)*Upcca_inv+(cos(T-2/3*pi)*COS-sin(T-2/3*pi)*SIN)*Upccb_inv+(cos(T+2/3*pi)*COS-sin(T+2/3*pi)*SIN)*Upccc_inv);
    uq_P_pu=-2/3*((sin(T)*COS+cos(T)*SIN)*Upcca_inv+(sin(T-2/3*pi)*COS+cos(T-2/3*pi)*SIN)*Upccb_inv+(sin(T+2/3*pi)*COS+cos(T+2/3*pi)*SIN)*Upccc_inv);
    #idq_N_pu
    id_N_pu= 2/3*((cos(-T)*COS_N-sin(-T)*SIN_N)*y36+(cos(-T-2/3*pi)*COS_N-sin(-T-2/3*pi)*SIN_N)*y37+(cos(-T+2/3*pi)*COS_N-sin(-T+2/3*pi)*SIN_N)*(-y36-y37));
    iq_N_pu=-2/3*((sin(-T)*COS_N+cos(-T)*SIN_N)*y36+(sin(-T-2/3*pi)*COS_N+cos(-T-2/3*pi)*SIN_N)*y37+(sin(-T+2/3*pi)*COS_N+cos(-T+2/3*pi)*SIN_N)*(-y36-y37));
    #udq_N_pu
    ud_N_pu= 2/3*((cos(-T)*COS_N-sin(-T)*SIN_N)*Upcca_inv+(cos(-T-2/3*pi)*COS_N-sin(-T-2/3*pi)*SIN_N)*Upccb_inv+(cos(-T+2/3*pi)*COS_N-sin(-T+2/3*pi)*SIN_N)*Upccc_inv);
    uq_N_pu=-2/3*((sin(-T)*COS_N+cos(-T)*SIN_N)*Upcca_inv+(sin(-T-2/3*pi)*COS_N+cos(-T-2/3*pi)*SIN_N)*Upccb_inv+(sin(-T+2/3*pi)*COS_N+cos(-T+2/3*pi)*SIN_N)*Upccc_inv);
    #idq_com
    id_com= 2/3*((cos(-2*T)*COS_2N-sin(-2*T)*SIN_2N)*y38+(cos(-2*T-2/3*pi)*COS_2N-sin(-2*T-2/3*pi)*SIN_2N)*y39+(cos(-2*T+2/3*pi)*COS_2N-sin(-2*T+2/3*pi)*SIN_2N)*y40);
    iq_com=-2/3*((sin(-2*T)*COS_2N+cos(-2*T)*SIN_2N)*y38+(sin(-2*T-2/3*pi)*COS_2N+cos(-2*T-2/3*pi)*SIN_2N)*y39+(sin(-2*T+2/3*pi)*COS_2N+cos(-2*T+2/3*pi)*SIN_2N)*y40);
    ##
    Tddsrf=0.004501581580786;
    dydt[49]=(1/Tddsrf)*(-y49+(ud_P_pu-( (cos(2*T)*COS_2-sin(2*T)*SIN_2)*y51+(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y52)));
    dydt[50]=(1/Tddsrf)*(-y50+(uq_P_pu-(-(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y51+(cos(2*T)*COS_2-sin(2*T)*SIN_2)*y52)));
    dydt[51]=(1/Tddsrf)*(-y51+(ud_N_pu-( (cos(-2*T)*COS_2N-sin(-2*T)*SIN_2N)*y49+(sin(-2*T)*COS_2N+cos(-2*T)*SIN_2N)*y50)));
    dydt[52]=(1/Tddsrf)*(-y52+(uq_N_pu-(-(sin(-2*T)*COS_2N+cos(-2*T)*SIN_2N)*y49+(cos(-2*T)*COS_2N-sin(-2*T)*SIN_2N)*y50)));
    ##
    dydt[53]=(1/Tddsrf)*(-y53+(id_P_pu-( (cos(2*T)*COS_2-sin(2*T)*SIN_2)*y55+(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y56)));
    dydt[54]=(1/Tddsrf)*(-y54+(iq_P_pu-(-(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y55+(cos(2*T)*COS_2-sin(2*T)*SIN_2)*y56)));
    dydt[55]=(1/Tddsrf)*(-y55+(id_N_pu-( (cos(-2*T)*COS_2N-sin(-2*T)*SIN_2N)*y53+(sin(-2*T)*COS_2N+cos(-2*T)*SIN_2N)*y54)));
    dydt[56]=(1/Tddsrf)*(-y56+(iq_N_pu-(-(sin(-2*T)*COS_2N+cos(-2*T)*SIN_2N)*y53+(cos(-2*T)*COS_2N-sin(-2*T)*SIN_2N)*y54)));
    #Phase Lock Loop
    # Kp_PLL_inv=10;
    # Ki_PLL_inv=1/0.02;
    dydt[57]=Ki_PLL_inv*(uq_P_pu-(-(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y51+(cos(2*T)*COS_2-sin(2*T)*SIN_2)*y52));
    dydt[58]=y57+Kp_PLL_inv*(uq_P_pu-(-(sin(2*T)*COS_2+cos(2*T)*SIN_2)*y51+(cos(2*T)*COS_2-sin(2*T)*SIN_2)*y52));
#outter loop 1:constant voltage
    Tm=0.0012;
    Gm=1/400;
    dydt[59]=-1/Tm*y59+Gm/Tm*(Ubase2_P2P_inv*udc_inv);
    # dydt[60]=Ki_udc*(Uref-y59);
    Idref_ori=Kp_udc*(Uref-y59)+y60;
    # Idref=idrefmax*((1+tanh(AA_ctrl*(Idref_ori-idrefmax)))/2)+idrefmin*((1+tanh(AA_ctrl*(idrefmin-Idref_ori)))/2)+(1-(((1+tanh(AA_ctrl*(Idref_ori-idrefmax)))/2)+((1+tanh(AA_ctrl*(idrefmin-Idref_ori)))/2)))*Idref_ori;
    
    
    ind_Id_Up=((1+tanh(AA_ctrl*(Idref_ori-idrefmax)))/2);
    ind_Id_Lo=((1+tanh(AA_ctrl*(idrefmin-Idref_ori)))/2);

    ind_Id_Up_inn=((1+tanh(AA_ctrl*(y60-(idrefmax-del))))/2);
    ind_Id_Lo_inn=((1+tanh(AA_ctrl*((idrefmin+del)-y60)))/2);

    Lim_I_id=-10*(y60-idrefmax)*ind_Id_Up_inn-10*(y60-idrefmin)*ind_Id_Lo_inn+(1-(ind_Id_Up_inn+ind_Id_Lo_inn))*Ki_udc*(Uref-y59);
    Ki_udc_rea=(1-(ind_Id_Up+ind_Id_Lo))*Ki_udc+(ind_Id_Up+ind_Id_Lo)*cof
    # Ki_udc_rea=Ki_udc
    # dydt[60]=(1-(ind_Id_Up+ind_Id_Lo))*Ki_udc*(Uref-y59)+(ind_Id_Up+ind_Id_Lo)*(-cof*y60);
    dydt[60]=Ki_udc_rea*(Uref-y59);
    Idref=idrefmax*ind_Id_Up+idrefmin*ind_Id_Lo+(1-(ind_Id_Up+ind_Id_Lo))*Idref_ori;
#outter loop 2:constant Q
    Q_ref= 0.2;
    # dydt[61]=Ki_cQ*(Q_ref-y49*y54);
    Iqref_ori=Kp_cQ*(Q_ref-y49*y54)+y61;
    # Iqref=1.5*((1+tanh(AA_ctrl*(Iqref_ori-1.5)))/2)-1.5*((1+tanh(AA_ctrl*(-1.5-Iqref_ori)))/2)+(1-(((1+tanh(AA_ctrl*(Iqref_ori-1.5)))/2)+((1+tanh(AA_ctrl*(-1.5-Iqref_ori)))/2)))*Iqref_ori;
    
    ind_Iq_Up=((1+tanh(AA_ctrl*(Iqref_ori-iqrefmax)))/2);
    ind_Iq_Lo=((1+tanh(AA_ctrl*(iqrefmin-Iqref_ori)))/2);

    ind_Iq_Up_inn=((1+tanh(AA_ctrl*(y61-(iqrefmax-del))))/2);
    ind_Iq_Lo_inn=((1+tanh(AA_ctrl*((iqrefmin+del)-y61)))/2);

    Lim_I_iq=-10*(y61-iqrefmax)*ind_Iq_Up_inn-10*(y61-iqrefmin)*ind_Iq_Lo_inn+(1-(ind_Iq_Up_inn+ind_Iq_Lo_inn))*Ki_cQ*(Q_ref-y49*y54);
    Ki_cQ_rea=(1-(ind_Iq_Up+ind_Iq_Lo))*Ki_cQ+(ind_Iq_Up+ind_Iq_Lo)*cof
    # Ki_cQ_rea=Ki_cQ
    # dydt[61]=(1-(ind_Iq_Up+ind_Iq_Lo))*Ki_cQ*(Q_ref-y49*y54)+(ind_Iq_Up+ind_Iq_Lo)*(-cof*y61);
    dydt[61]=Ki_cQ_rea*(Q_ref-y49*y54);
    Iqref=iqrefmax*ind_Iq_Up+iqrefmin*ind_Iq_Lo+(1-(ind_Iq_Up+ind_Iq_Lo))*Iqref_ori;
#Positive inner loop
    dydt[62]=Ki_CC*(Idref-y53);
    dydt[63]=Ki_CC*(Iqref-y54);
    dydt[64]=Ki_CC*(0-y55);
    dydt[65]=Ki_CC*(0-y56);
    
    Udref_P=y49-(y62+Kp_CC*(Idref-y53))+Kd*y54;
    Uqref_P=y50-(y63+Kp_CC*(Iqref-y54))-Kd*y53;
    Udref_N=y51-(y64+Kp_CC*(0-y55))-Kd*y56;
    Uqref_N=y52-(y65+Kp_CC*(0-y56))+Kd*y55;
#Circul
    Kp_CCSC=20;Ki_CCSC=1/0.01;
    dydt[66]=Ki_CCSC*(0-Ibase2_P2P_inv*id_com);
    dydt[67]=Ki_CCSC*(0-Ibase2_P2P_inv*iq_com);
    Udref_CCSC=Kp_CCSC*(0-Ibase2_P2P_inv*id_com)+y66;
    Uqref_CCSC=Kp_CCSC*(0-Ibase2_P2P_inv*iq_com)+y67;
#
    Uaref=(COS*cos(T)-SIN*sin(T))*Udref_P+(-(SIN*cos(T)+COS*sin(T)))*Uqref_P+(COS_N*cos(-T)-SIN_N*sin(-T))*Udref_N+(-(SIN_N*cos(-T)+COS_N*sin(-T)))*Uqref_N;
    Ubref=(COS*cos(T-2*pi/3)-SIN*sin(T-2*pi/3))*Udref_P+(-(SIN*cos(T-2*pi/3)+COS*sin(T-2*pi/3)))*Uqref_P+(COS_N*cos(-T-2*pi/3)-SIN_N*sin(-T-2*pi/3))*Udref_N+(-(SIN_N*cos(-T-2*pi/3)+COS_N*sin(-T-2*pi/3)))*Uqref_N;
    Ucref=(COS*cos(T+2*pi/3)-SIN*sin(T+2*pi/3))*Udref_P+(-(SIN*cos(T+2*pi/3)+COS*sin(T+2*pi/3)))*Uqref_P+(COS_N*cos(-T+2*pi/3)-SIN_N*sin(-T+2*pi/3))*Udref_N+(-(SIN_N*cos(-T+2*pi/3)+COS_N*sin(-T+2*pi/3)))*Uqref_N;
    
    Uaref_CCSC=(COS_2N*cos(-2*T)-SIN_2N*sin(-2*T))*Udref_CCSC+(-(SIN_2N*cos(-2*T)+COS_2N*sin(-2*T)))*Uqref_CCSC;
    Ubref_CCSC=(COS_2N*cos(-2*T-2*pi/3)-SIN_2N*sin(-2*T-2*pi/3))*Udref_CCSC+(-(SIN_2N*cos(-2*T-2*pi/3)+COS_2N*sin(-2*T-2*pi/3)))*Uqref_CCSC;
    Ucref_CCSC=(COS_2N*cos(-2*T+2*pi/3)-SIN_2N*sin(-2*T+2*pi/3))*Udref_CCSC+(-(SIN_2N*cos(-2*T+2*pi/3)+COS_2N*sin(-2*T+2*pi/3)))*Uqref_CCSC;
    
    
    na_T_inter=(0.5*400-Ubase2_P2P_inv*Uaref-Uaref_CCSC)/2.1/200;
    na_B_inter=(0.5*400+Ubase2_P2P_inv*Uaref-Uaref_CCSC)/2.1/200;
    nb_T_inter=(0.5*400-Ubase2_P2P_inv*Ubref-Ubref_CCSC)/2.1/200;
    nb_B_inter=(0.5*400+Ubase2_P2P_inv*Ubref-Ubref_CCSC)/2.1/200;
    nc_T_inter=(0.5*400-Ubase2_P2P_inv*Ucref-Ucref_CCSC)/2.1/200;
    nc_B_inter=(0.5*400+Ubase2_P2P_inv*Ucref-Ucref_CCSC)/2.1/200;
    
    na_T=na_T_inter;
    # ((1+tanh(AA_ctrl*(na_T_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-na_T_inter)))/2)+(1-(((1+tanh(AA_ctrl*(na_T_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-na_T_inter)))/2)))*na_T_inter;
    na_B=na_B_inter;
    # ((1+tanh(AA_ctrl*(na_B_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-na_B_inter)))/2)+(1-(((1+tanh(AA_ctrl*(na_B_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-na_B_inter)))/2)))*na_B_inter;
    nb_T=nb_T_inter;
    # ((1+tanh(AA_ctrl*(nb_T_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-nb_T_inter)))/2)+(1-(((1+tanh(AA_ctrl*(nb_T_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-nb_T_inter)))/2)))*nb_T_inter;
    nb_B=nb_B_inter;
    # ((1+tanh(AA_ctrl*(nb_B_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-nb_B_inter)))/2)+(1-(((1+tanh(AA_ctrl*(nb_B_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-nb_B_inter)))/2)))*nb_B_inter;
    nc_T=nc_T_inter;
    # ((1+tanh(AA_ctrl*(nc_T_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-nc_T_inter)))/2)+(1-(((1+tanh(AA_ctrl*(nc_T_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-nc_T_inter)))/2)))*nc_T_inter;
    nc_B=nc_B_inter;
    # ((1+tanh(AA_ctrl*(nc_B_inter-1)))/2)-0.7*((1+tanh(AA_ctrl*(-0.7-nc_B_inter)))/2)+(1-(((1+tanh(AA_ctrl*(nc_B_inter-1)))/2)+((1+tanh(AA_ctrl*(-0.7-nc_B_inter)))/2)))*nc_B_inter;
    ##main circuit
    ua_tri=(na_B*y42-na_T*y41)/2;ub_tri=(nb_B*y44-nb_T*y43)/2;uc_tri=(nc_B*y46-nc_T*y45)/2;
    unN=(1/3)*((Upcca_inv+Upccb_inv+Upccc_inv)-(ua_tri+ub_tri+uc_tri));
    # unN=(1/3)*(-(ua_tri+ub_tri+uc_tri));
    dydt[36]=(1/Leq)*(-unN+Upcca_inv-ua_tri);
    dydt[37]=(1/Leq)*(-unN+Upccb_inv-ub_tri);
    ####
    ##Common mode Current
    dydt[38]=((1/Larm)*(udc_inv/2-((na_B*y42+na_T*y41)/2)));
    dydt[39]=((1/Larm)*(udc_inv/2-((nb_B*y44+nb_T*y43)/2)));
    dydt[40]=((1/Larm)*(udc_inv/2-((nc_B*y46+nc_T*y45)/2)));
    ##Sum of sm voltage
    dydt[41]= (1/Ceq)*na_T*((2*(y38)-((y36)))/2);
    dydt[42]= (1/Ceq)*na_B*((2*(y38)+((y36)))/2);
    dydt[43]= (1/Ceq)*nb_T*((2*(y39)-((y37)))/2);
    dydt[44]= (1/Ceq)*nb_B*((2*(y39)+((y37)))/2);
    dydt[45]=(1/Ceq)*nc_T*((2*(y40)-(-(y36)-(y37)))/2);
    dydt[46]=(1/Ceq)*nc_B*((2*(y40)+(-(y36)-(y37)))/2);
##
    dydt[47]=(Rg_inv/Lg_inv)*((y36+y73)-y47);
    dydt[48]=(Rg_inv/Lg_inv)*(y37-y48);
    dydt[72]=(Rg_inv/Lg_inv)*((-y36-y37)-y72);

    # dydt[47]=(Rg_inv/Lg_inv)*(y36-y47);
    # dydt[48]=(Rg_inv/Lg_inv)*(y37-y48);

##
    # Lfault_inv=Lfault_real*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    # Lfault_inv=100*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    # Rfault_inv=1e-3*Tra_inv/(Ubase2_L2L_inv^2/Sbase_inv);
    # dydt[73]=1/Lfault_inv*(Upcca_inv);
    dydt[73]=-10*y73;

##
    dydt[68]=(1/Cdc_inv)*(-y38-y39-y40-y69);
    dydt[69]=(1/Ldc_inv)*(y68-Edc_inv);
    # ## 
    dydt[70]=1/Cline*(idc_rec_end+idc_inv_end);
##
    dydt
    end
2

Are the stats significantly different? Also, what versions of OrdinaryDiffEq OrdinaryDiffEqRosenbrock, and OrdinaryDiffEqDifferentiation? that is a lot more important for figuring out what’s going wrong.

3

7.14.0：
name = “OrdinaryDiffEqDifferentiation”
version = “1.3.0”

name = “OrdinaryDiffEqRosenbrock”
version = “1.4.0”

name = “OrdinaryDiffEq”
version = “6.90.1”

7.9.0：I only find OrdinaryDiffEq in C:.….julia\packages
name = “OrdinaryDiffEq”
version = “6.58.2”

1 Like
4

You’re testing with an auto-switch method. Test with just the Vern6 and just the Rodas5P and see which one it’s matching.

5

the other possibility is that there method chosen changed. I switched the stiffness calculation on the nonstiff side from a 2 norm to an inf norm and it’s possible that doing so makes it use a different algorithm

6

Hi, Oscar and Chris, Here are the test results.
DifferentialEquation@7.9.0:

AutoVern6(Rodas5P(),nonstifftol = 15/ 10):
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.418 s …    2.346 s  ┊ GC (min … max):  0.00% … 41.36%
Memory estimate: 661.43 MiB, allocs estimate: 2418966.

Vern6():
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.451 s …    1.982 s  ┊ GC (min … max):  0.00% … 26.63%
Memory estimate: 659.45 MiB, allocs estimate: 2418817.

Rodas5P():
BenchmarkTools.Trial: 4 samples with 1 evaluation.
Range (min … max):  22.905 s …   24.300 s  ┊ GC (min … max): 2.90% … 2.82%
Memory estimate: 1.88 GiB, allocs estimate: 21146606.

DifferentialEquation@7.14.0:

AutoVern6(Rodas5P(),nonstifftol = 15/ 10):
BenchmarkTools.Trial: 2 samples with 1 evaluation.
Range (min … max):  20.355 s …  20.452 s  ┊ GC (min … max): 7.00% … 7.22%
 Memory estimate: 1.52 GiB, allocs estimate: 15735599.

Vern6():
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.651 s …    2.611 s  ┊ GC (min … max):  9.16% … 44.33%
Memory estimate: 660.53 MiB, allocs estimate: 3311381.

Rodas5P():
BenchmarkTools.Trial: 1 sample with 1 evaluation.
Single result which took 26.859 s (5.53% GC) to evaluate,
with a memory estimate of 1.81 GiB, over 19996952 allocations.

It seems that 7.9.0 automatically picks Vern6, while 7.14.0 goes with Rodas5P, which is pretty strange. Another thing I’m really curious about is why Vern6 for non-stiff problems is more efficient. The plot below shows the waveforms of this system in steady state, and the cornered lines for y19 and y20 actually represent state-dependent switching, which should make the problem pretty stiff. Maybe I can stick with Vern6, but are there any solvers that might work better for this kind of issue?"

untitled1
untitled1840×630 22.6 KB

untitled2
untitled2840×630 17.2 KB

7

That would be norm change. It might’ve become a bit too sensitive. Can you open an issue on OrdinaryDiffEq.jl?

8

Sure! But are there any algorithms or parameter settings that can solve this problem more efficiently?

9

You didn’t show any events? And this doesn’t look stiff at all? From the looks of it, Vern6 or 7 seem like a good fit. Otherwise try VCABM?

10

I’m not quite sure what you mean by “event.” If you’re referring to discrete events introduced by switching (which can be represented using the Heaviside function), they’ve already been approximated into a continuous form. Specifically, Heaviside(x) is approximated as h(x)=0.5×(1+tanh⁡(kx))·.
You can check out the link below for more details:

From an intuitive perspective, the solver may spends a lot of time around regions where x is near zero (similar to those corners in the graph for y19 and y20), which could make the problem stiff? Anyway, the Vern algorithm looks pretty good, and VCABM seems to be even more efficient. Thank you so much for your explanation!

1 Like
11

Yeah that’s not very likely to be that stiff. You’d need eigenvalue differences of around >1e6 for an implicit method to make sense, the overhead of the methods that add stability is fairly high and so if you have this kind of “mild stiffness” then it usually just makes sense to take an explicit method and just push it right through with a few extra steps. Or something in the middle, like a ROCK method (ROCK2) or VCABM, which is an implicit method with less stable implicit handling and adaptive order to the stability region, would likely do better in most cases that look like this.