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

**URL:** <https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377>\
**Category:** Numerics\
**Tags:** package, differentialequation\
**Created:** [January 2, 2025, 4:38pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377 "2025-01-02T16:38:52Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![equi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/equi/32/214471_2.png) [@equi](https://discourse.julialang.org/u/equi)\
**Post date:** [January 2, 2025, 4:38pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/1 "2025-01-02T16:38:52Z")

</div>

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:

```julia
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:

```julia
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
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:

```julia
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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 2, 2025, 5:31pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/2 "2025-01-02T17:31:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![equi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/equi/32/214471_2.png) [@equi](https://discourse.julialang.org/u/equi)\
**Post date:** [January 2, 2025, 6:34pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/3 "2025-01-02T18:34:04Z")

</div>

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”

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 3, 2025, 1:03am UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/4 "2025-01-03T01:03:19Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 3, 2025, 1:57am UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/5 "2025-01-03T01:57:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![equi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/equi/32/214471_2.png) [@equi](https://discourse.julialang.org/u/equi)\
**Post date:** [January 3, 2025, 3:42pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/6 "2025-01-03T15:42:36Z")

</div>

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

```julia
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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/4/6/463fe98303437bfdc4258bab732bd85b77c1d830.jpeg)  
 ![untitled2](https://global.discourse-cdn.com/julialang/original/3X/2/2/227bb80de993c4ce74ff6f6b6e5ef80789122230.jpeg)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 4, 2025, 3:03am UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/7 "2025-01-04T03:03:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![equi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/equi/32/214471_2.png) [@equi](https://discourse.julialang.org/u/equi)\
**Post date:** [January 4, 2025, 2:34pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/8 "2025-01-04T14:34:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 4, 2025, 4:39pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/9 "2025-01-04T16:39:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![equi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/equi/32/214471_2.png) [@equi](https://discourse.julialang.org/u/equi)\
**Post date:** [January 5, 2025, 2:18pm UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/10 "2025-01-05T14:18:36Z")

</div>

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:

> **[Approximating piecewise nonlinearities in dynamic systems with sigmoid...](https://link.springer.com/article/10.1007/s11071-023-08293-1)**
>
> In the industry field, the increasingly stringent requirements of lightweight structures are exposing the ultimately nonlinear nature of mechanical systems. This is extremely true for systems with moving parts and loose fixtures which show piecewise...

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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 6, 2025, 10:12am UTC](https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/11 "2025-01-06T10:12:14Z")

</div>

> [@equi](#):
>
> 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))`·.

> [@equi](#):
>
> 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?

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.