Solve DAEs using IDASolve function of SUNDIALS 2.6.2 version

I am trying to solve the system of differential algebraic equations which are used to model the IEEE 9 bus system. The differential equations determines the dynamic behaviour of the generators while the algebraic equations consists of the stator and the network (power flow) equations.

I have attached a MATLAB code in which the function "fun_DAE" gives out all of the equations in the form of a column vectors. ' F_D ' gives the 21 differential equations (7 for each generator), ' F_SA ' gives 18 network equations (real and reactive power balances for 9 buses) and ' F_SA' gives the 6 stator equations (2 for each generator). ' F_SA' and ' F_N ' constitutes the algebraic constraints.

the differential variables are for each generator. The initial values of this state variables are calculated and stored as matrix.

I need some help in using the IDAInit, IDASetOptions and IDASolve functions or any other relevant functions to be used to get a solution.

Thank You.

[Torsten]

Could you include your code as an .m-file so that we are able to run it here directly ?

[Ajinkya]

@Torsten

" fun_DAE.m" is function containing the DAEs while "YBus" gives the Ybus matrix and solves the load flow giving "V" voltages at the buses. The reference power system chosen is the IEEE 9 bus system.

• "WSCC_9_Bus.xlsx" constains the system data
• first need to run the "YBus()" code to get the [Y, V]
• get inside the "MainCode.m" and run it to get the values that are needed for "fun_DAE". (inside the .m file I have tried solving the system using the Newtons method but it was a fail.

[Torsten]

"fun_DAE" returns a 45x1 column vector. I don't know what the 45 unknowns are in your input list to "fun_DAE" that correspond to these 45 equations.

Maybe you have a document in which you list the 45 equations together with the 45 corresponding unknowns.

[Ajinkya]

In the attached pdf are the actual equations. There are only 21 differential equations while rest 26 are the algebraic equations.

[Torsten]

What are the unknowns in the algebraic equations ?

[Ajinkya]

Id1 Id2 Id3 Iq1 Iq2 Iq3 are the unknowns.

The voltages V1 to V9 and its angles theta1 to theta9 are all known values

[Torsten]

And what are the last 18 algebraic equations for that you list in your "Dae.pdf" ? It seems they are not needed in your solution process then.

[Ajinkya]

Those 18 equations are the power balance equations at each bus.

First 9 equations are Real power equations while the rest 9 are reactive power equations.

I have modified the fun_DAE function to have three arguments.

I am setting the ode and assigning the initial values using vector "x" which is a (21x1) vector.

F = ode;
F.ODEFcn = @fun_DAE;

Do I need to exclusively pass the " x " vector into the function so that the ode object recognizes it as the solution. Because there are many variables considered inside the function as parameters of the DAE. I am unable to proceed further.

% dataM, dataE are the structures defined to contain the Machine & Exciter
% data
% Init structure contains differential variables initial values and other
% data
function F = fun_DAE(dataM,dataE,Init)

% MACHINE DATA    EXTRACTING FROM THE STRUCT 
        [DM, H, M, Xd, Xdd] = deal(dataM.DM, dataM.H, dataM.M, dataM.Xd,dataM.Xdd);
     [Xq, Xqq, Td0, Tq0] = deal(dataM.Xq,dataM.Xqq, dataM.Td0, dataM.Tq0);
        D = DM.*M;
    ws = 2*pi*60;
% EXCITER DATA    EXTRACTING FROM THE STRUCT 
    [KA, TA, KE, TE, KF, TF] = deal(dataE.KA, dataE.TA, dataE.KE, dataE.TE, dataE.KF, dataE.TF);
%  x IS THE DIFFERENTIAL VARIABLE'S INITIAL VALUES VECTOR (21x1)
    x = Init.x;
    Id = Init.Idq(1,:);
    Iq = Init.Idq(2,:);
    Vref = Init.Vref_Tm(1,:);
    Tm = Init.Vref_Tm(2,:);
    V = Init.V;
    th = Init.th;
    Y = Init.Y;
    % WRITING DIFFERENTIAL ALGEBRAIC EQUATIONS
% GENERATOR 1
    F11 = -x(1,1)/Td0(1) - (Xd(1)-Xdd(1))*Id(1)/Td0(1) + x(5,1);
    F12 = -x(2,1)/Td0(1) - (Xq(1)-Xqq(1))*Iq(1);
    F13 = x(4,1) - ws;
    F14 = (Tm(1) - x(2,1)*Id(1) - x(1,1)*Iq(1) -(Xqq(1)-Xdd(1))*Id(1)*Iq(1))/M(1) - D(1)*(x(4,1) -ws);
    F15 = -(KE(1) + 0.0039*exp(1.555*x(5,1)))*x(5,1)/TE(1) + x(6,1)/TE(1);
    F16 = -x(7,1)/TF(1) + KF(1)*x(5,1)/((TF(1))^2);
    F17 = (-x(6,1) + KA(1)*x(7,1) - (KA(1)*KF(1)*x(5,1))/TF(1) + KA(1)*(V(1)-Vref(1)))/TA(1);
% GENERATOR 2
    F21 = -x(1,2)/Td0(2) - (Xd(2)-Xdd(2))*Id(2)/Td0(2) + x(5,2);
    F22 = -x(2,2)/Td0(2) - (Xq(2)-Xqq(2))*Iq(2);
    F23 = x(4,2) - ws;
    F24 = (Tm(2) - x(2,2)*Id(2) - x(1,2)*Iq(2) -(Xqq(2)-Xdd(2))*Id(2)*Iq(2))/M(2) - D(2)*(x(4,2) -ws);
    F25 = -(KE(2) + 0.0039*exp(1.555*x(5,2)))*x(5,2)/TE(2) + x(6,2)/TE(2);
    F26 = -x(7,2)/TF(2) + KF(2)*x(5,2)/((TF(2))^2);
    F27 = (-x(6,2) + KA(2)*x(7,2) - (KA(2)*KF(2)*x(5,2))/TF(2) + KA(2)*(V(2)-Vref(2)))/TA(2);
% GENERATOR 3
    F31 = -x(1,3)/Td0(3) - (Xd(3)-Xdd(3))*Id(3)/Td0(3) + x(5,3);
    F32 = -x(2,3)/Td0(3) - (Xq(3)-Xqq(3))*Iq(3);
    F33 = x(4,3) - ws;
    F34 = (Tm(3) - x(2,3)*Id(3) - x(1,3)*Iq(3) -(Xqq(3)-Xdd(3))*Id(3)*Iq(3))/M(3) - D(3)*(x(4,3) -ws);
    F35 = -(KE(3) + 0.0039*exp(1.555*x(5,3)))*x(5,3)/TE(3) + x(6,3)/TE(3);
    F36 = -x(7,3)/TF(3) + KF(3)*x(5,3)/((TF(3))^2);
    F37 = (-x(6,3) + KA(3)*x(7,3) - (KA(3)*KF(3)*x(5,3))/TF(3) + KA(3)*(V(3)-Vref(3)))/TA(3);
% CONSOLIDATING ALL EQUATIONS
    F_D = [F11;F12;F13;F14;F15;F16;F17;F21;F22;F23;F24;F25;F26;F27;F31;F32;F33;F34;F35;F36;F37];

% STATOR ALGEBRAIC EQUATIONS
    F18 = Id(1)*Xd(1) + V(1)*cos(x(3,1)-th(1)) - x(1,1);
    F19 = Iq(1)*Xq(1) - V(1)*sin(x(3,1)-th(1)) - x(2,1);
    F28 = Id(2)*Xd(2) + V(2)*cos(x(3,2)-th(2)) - x(1,2);
    F29 = Iq(2)*Xq(2) - V(2)*sin(x(3,2)-th(2)) - x(2,2);
    F38 = Id(3)*Xd(3) + V(3)*cos(x(3,3)-th(3)) - x(1,3);
    F39 = Iq(3)*Xq(3) - V(3)*sin(x(3,3)-th(3)) - x(2,3);
    
    F_SA = [F18;F19;F28;F29;F38;F39];

% NETWORK ALGEBRAIC EQUATIONS 
% REAL POWER, 9 EQUATIONS FOR 9 BUSES
    F110 = Id(1)*V(1)*sin(x(3,1)-th(1)) + Iq(1)*V(1)*cos(x(3,1)-th(1)) - 17.36*V(1)*V(4)*sin(th(1)-th(4));
    F210 = Id(2)*V(2)*sin(x(3,2)-th(2)) + Iq(2)*V(2)*cos(x(3,2)-th(2)) - 16*V(2)*V(7)*sin(th(2)-th(7));
    F310 = Id(3)*V(3)*sin(x(3,3)-th(3)) + Iq(3)*V(3)*cos(x(3,3)-th(3)) - 17.06*V(3)*V(9)*sin(th(3)-th(9));
    F410 =  -abs(Y(4,1))*V(1)*V(4)*sin(th(4)-th(1)) - abs(Y(4,4))*V(4)^2*cos(angle(Y(4,4))) -abs(Y(4,5))*V(4)*V(5)*cos(th(4)-th(5)-angle(Y(4,5))) -abs(Y(4,6))*V(4)*V(6)*cos(th(4)-th(6)-angle(Y(4,6)));
    F510 =  -1.25 - abs(Y(5,4))*V(5)*V(4)*cos(th(5)-th(4)-angle(Y(5,4))) - abs(Y(5,5))*V(5)^2*cos(angle(Y(5,5))) -abs(Y(5,7))*V(5)*V(7)*cos(th(5)-th(7)-angle(Y(5,7)));
    F610 =  -0.9 - abs(Y(6,4))*V(6)*V(4)*cos(th(6)-th(4)-angle(Y(6,4))) - abs(Y(6,6))*V(6)^2*cos(angle(Y(6,6))) -abs(Y(6,9))*V(6)*V(9)*cos(th(6)-th(9)-angle(Y(6,9)));
    F710 =  -abs(Y(7,2))*V(7)*V(2)*sin(th(7)-th(2)) - abs(Y(7,7))*V(7)^2*cos(angle(Y(7,7))) -abs(Y(7,5))*V(7)*V(5)*cos(th(7)-th(5)-angle(Y(7,5))) -abs(Y(7,8))*V(7)*V(8)*cos(th(7)-th(8)-angle(Y(7,8)));
    F810 =  -1 -abs(Y(8,7))*V(7)*V(8)*cos(th(8)-th(7)-angle(Y(8,7))) - abs(Y(8,8))*V(8)^2*cos(angle(Y(8,8))) -abs(Y(8,9))*V(9)*V(8)*cos(th(8)-th( 9)-angle(Y(8,9)));
    F910 =  -abs(Y(9,3))*V(9)*V(3)*sin(th(9)-th(3)) - abs(Y(9,9))*V(9)^2*cos(angle(Y(9,9))) -abs(Y(9,6))*V(9)*V(6)*cos(th(9)-th(6)-angle(Y(9,6))) -abs(Y(9,8))*V(9)*V(8)*cos(th(9)-th(8)-angle(Y(9,8)));
% RECTIVE POWER
    F111 = Id(1)*V(1)*cos(x(3,1)-th(1)) - Iq(1)*V(1)*sin(x(3,1)-th(1)) + 17.36*V(1)*V(4)*cos(th(1)-th(4)) - 17.36*(V(1))^2;
    F211 = Id(2)*V(2)*cos(x(3,2)-th(2)) - Iq(2)*V(2)*sin(x(3,2)-th(2)) + 16*V(2)*V(7)*cos(th(2)-th(7)) - 16*(V(2))^2;
    F311 = Id(3)*V(3)*cos(x(3,3)-th(3)) - Iq(3)*V(3)*sin(x(3,3)-th(3)) + 17.06*V(3)*V(9)*cos(th(3)-th(9)) - 16*(V(3))^2;
    F411 = abs(Y(4,1))*V(1)*V(4)*cos(th(4)-th(1)) + abs(Y(4,4))*V(4)^2*sin(angle(Y(4,4))) -abs(Y(4,5))*V(4)*V(5)*sin(th(4)-th(5)-angle(Y(4,5))) -abs(Y(4,6))*V(4)*V(6)*sin(th(4)-th(6)-angle(Y(4,6)));
    F511 = -0.5 - abs(Y(5,4))*V(5)*V(4)*sin(th(5)-th(4)-angle(Y(5,4))) + abs(Y(5,5))*V(5)^2*sin(angle(Y(5,5))) -abs(Y(5,7))*V(5)*V(7)*sin(th(5)-th(7)-angle(Y(5,7)));
    F611 = -0.3 - abs(Y(6,4))*V(6)*V(4)*sin(th(6)-th(4)-angle(Y(6,4))) + abs(Y(6,6))*V(6)^2*sin(angle(Y(6,6))) -abs(Y(6,9))*V(6)*V(9)*sin(th(6)-th(9)-angle(Y(6,9)));
    F711 = abs(Y(7,2))*V(7)*V(2)*cos((th(7)-th(2)) + abs(Y(7,7))*V(7)^2*sin(angle(Y(7,7))) -abs(Y(7,5))*V(7)*V(5)*sin( (th(7)-th(5)-angle(Y(7,5)))) -abs(Y(7,8))*V(7)*V(8)*sin(th(7)-th(8)-angle(Y(7,8))));
    F811 = -0.35 -abs(Y(8,7))*V(7)*V(8)*sin(th(8)-th(7)-angle(Y(8,7))) + abs(Y(8,8))*V(8)^2*sin(angle(Y(8,8))) -abs(Y(8,9))*V(9)*V(8)*sin(th(8)-th(9)-angle(Y(8,9)));
    F911 = abs(Y(9,3))*V(9)*V(3)*cos(th(9)-th(3)) + abs(Y(9,9))*V(9)^2*sin(angle(Y(9,9))) -abs(Y(9,6))*V(9)*V(6)*sin(th(9)-th(6)-angle(Y(9,6))) -abs(Y(9,8))*V(9)*V(8)*sin(th(9)-th(8)-angle(Y(9,8)));

    F_N = [F110;F210;F310;F410;F510;F610;F710;F810;F910;F111;F211;F311;F411;F511;F611;F711;F811;F911];
    F = [F_D; F_SA; F_N];
end

[Ajinkya]

% this "data" MAT file contains the dataM & dataE struc
load data

% this "Bus" MAT file contains the "V" ,"th", "Y"
load Bus

% "Init" MAT file contains the initial values of Differential variables
% along with other parameters.
load Init

the following code is the main code that will initilize the DAE system. (Added this code for the reference only as

the above MAT files already contains the data that has been computed using the code below).

% MACHINE DATA
   [DM, H, M, Xd, Xdd, Xq, Xqq, Td0, Tq0] = deal(dataM.DM, dataM.H, dataM.M, dataM.Xd, ...
      dataM.Xdd, dataM.Xq,dataM.Xqq, dataM.Td0, dataM.Tq0);
   D = DM.*M;
   ws = 2*pi*60;
   % EXCITER DATA
    [KA, TA, KE, TE, KF, TF] = deal(dataE.KA, dataE.TA, dataE.KE, dataE.TE, dataE.KF, dataE.TF);
    
% DETERMINING THE INITIAL VALUES OF THE STATE VARIABLES
for k = 1:3
    Pg(k) = real(Sbus(k));
    Qg(k) = imag(Sbus(k));
    Ig(k) = conj((Pg(k) + 1i*Qg(k))/(V(k)*exp(1i*th(k))));
    delt(k) = angle(V(k)*exp(1i*th(k)) + 1i*Xq(k)*Ig(k));
    Idq(k) = abs(Ig(k))*exp(1i*(angle(Ig(k))-delt(k) + pi/2));
    Id(k) = real(Idq(k));   Iq(k) = imag(Idq(k));
    Vdq(k) = V(k)*exp(1i*(th(k) - delt(k) + pi/2));
    Ed(k) = real(Vdq(k) - Xqq(k)*imag(Idq(k)));
    Eq(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k));
    Efd(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k));
    Rf(k) = KF(k)*(imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k)))/TF(k);
    Vr(k) = (KE(k) + 0.0039*exp(1.555*Efd(k)))*Efd(k);
    Vref(k) = V(k) + Vr(k)/KA(k);
    Tm(k) = Ed(k)*real(Idq(k)) + Eq(k)*imag(Idq(k)) + (Xqq(k) - Xdd(k))*real(Idq(k))*imag(Idq(k));
    w(k) = ws;
end
for i = 1:3     % STATE VARIABLE VECTOR 
    %           1       2       3           4       5       6       7
    x(:,i) = [Eq(i); Ed(i); delt(i);    w(i);   Efd(i);     Vr(i); Rf(i)];
end
 x = [x(:,1);x(:,2);x(:,3)];
    Init.x = x;
    Init.Idq = [Id;Iq]; 
    Init.Vref_Tm = [Vref; Tm];
    Init.V = V;
    Init.th = th;
    Init.Y = Y;
    save Init Init 

[Torsten]

"fun_DAE" must have the form

function dydt = fun_DAE(t,y,data1,data2,....)

where data1, data2, ... can be any predefined parameters of the DAE system that you need to compute the time derivatives of the 21 differential equations.

The call to the ODE integrator should be

F = ode;
F.ODEFcn = @(t,y) fun_DAE(t,y,data1,data2,...);
F.InitialValue = [y1_0;y2_0;y3_0;...;y21_0];
F.Solver = 'idas';

where y1_0 is the initial value of E_q1, y2_0 is the initial value of E_d1, ..., y7_0 is the initial value of E_q2,... and so on.

Since I_q1,I_q2,I_q3,I_d1,I_d2,I_d3 can easily be computed from the stator equations, you don't need to define them as solution variables.

So as far as I see, fun_DAE only needs this part

function dydt = fun_DAE(t,y,data1,data2,...)
% GENERATOR 1
    F11 = -x(1,1)/Td0(1) - (Xd(1)-Xdd(1))*Id(1)/Td0(1) + x(5,1);
    F12 = -x(2,1)/Td0(1) - (Xq(1)-Xqq(1))*Iq(1);
    F13 = x(4,1) - ws;
    F14 = (Tm(1) - x(2,1)*Id(1) - x(1,1)*Iq(1) -(Xqq(1)-Xdd(1))*Id(1)*Iq(1))/M(1) - D(1)*(x(4,1) -ws);
    F15 = -(KE(1) + 0.0039*exp(1.555*x(5,1)))*x(5,1)/TE(1) + x(6,1)/TE(1);
    F16 = -x(7,1)/TF(1) + KF(1)*x(5,1)/((TF(1))^2);
    F17 = (-x(6,1) + KA(1)*x(7,1) - (KA(1)*KF(1)*x(5,1))/TF(1) + KA(1)*(V(1)-Vref(1)))/TA(1);
% GENERATOR 2
    F21 = -x(1,2)/Td0(2) - (Xd(2)-Xdd(2))*Id(2)/Td0(2) + x(5,2);
    F22 = -x(2,2)/Td0(2) - (Xq(2)-Xqq(2))*Iq(2);
    F23 = x(4,2) - ws;
    F24 = (Tm(2) - x(2,2)*Id(2) - x(1,2)*Iq(2) -(Xqq(2)-Xdd(2))*Id(2)*Iq(2))/M(2) - D(2)*(x(4,2) -ws);
    F25 = -(KE(2) + 0.0039*exp(1.555*x(5,2)))*x(5,2)/TE(2) + x(6,2)/TE(2);
    F26 = -x(7,2)/TF(2) + KF(2)*x(5,2)/((TF(2))^2);
    F27 = (-x(6,2) + KA(2)*x(7,2) - (KA(2)*KF(2)*x(5,2))/TF(2) + KA(2)*(V(2)-Vref(2)))/TA(2);
% GENERATOR 3
    F31 = -x(1,3)/Td0(3) - (Xd(3)-Xdd(3))*Id(3)/Td0(3) + x(5,3);
    F32 = -x(2,3)/Td0(3) - (Xq(3)-Xqq(3))*Iq(3);
    F33 = x(4,3) - ws;
    F34 = (Tm(3) - x(2,3)*Id(3) - x(1,3)*Iq(3) -(Xqq(3)-Xdd(3))*Id(3)*Iq(3))/M(3) - D(3)*(x(4,3) -ws);
    F35 = -(KE(3) + 0.0039*exp(1.555*x(5,3)))*x(5,3)/TE(3) + x(6,3)/TE(3);
    F36 = -x(7,3)/TF(3) + KF(3)*x(5,3)/((TF(3))^2);
    F37 = (-x(6,3) + KA(3)*x(7,3) - (KA(3)*KF(3)*x(5,3))/TF(3) + KA(3)*(V(3)-Vref(3)))/TA(3);
% CONSOLIDATING ALL EQUATIONS
    dydt = [F11;F12;F13;F14;F15;F16;F17;F21;F22;F23;F24;F25;F26;F27;F31;F32;F33;F34;F35;F36;F37];
end

where I don't understand which of the variables you use are known in advance and which are the unknowns y.

[Ajinkya]

The variable " x " is the unknown which will be computed.

Rest of the variables are either system constants or are already calculated in the code.

I will try running this code once.

Thank you so much for your time and help !!

[Torsten]

The initial conditions must be given in a vector of length 21, not in a matrix of size 7x3. The vector y you receive from "idas" in function "funDAE" is also of size 21x1, not 7x3. You should keep this in mind when reorganizing your code.

[Ajinkya]

You said that I don't need to include the last 18 equations but those will constitute the algebraic eqns in the DAE .

[Torsten]

This was your answer when I asked what the unknowns are in the algebraic equations:

Id1 Id2 Id3 Iq1 Iq2 Iq3 are the unknowns.

The voltages V1 to V9 and its angles theta1 to theta9 are all known values

This would mean that you only need to know Id1 Id2 Id3 Iq1 Iq2 Iq3. And these values can easily be computed directly in "funDAE" without defining them as unknowns for the ODE integrator:

Id1 = (Eq1 - V1*cos(delta1-theta1))/0.0608 
Id2 = (Eq2 - V2*cos(delta2-theta2))/0.1198     
Id3 = (Eq3 - V3*cos(delta3-theta3))/(-0.0608)    
Iq1 = (Ed1 - V1*cos(delta1-theta1))/(-0.0969)
Iq2 = (Ed2 - V2*cos(delta2-theta2))/(-0.1969)
Iq3 = (Ed3 - V3*cos(delta3-theta3))/(-0.0969)

[Ajinkya]

Oh sorry, my bad for this one !!

[Torsten]

So there are additional variables in the algebraic equations that are unknown ? Or how to interprete "Oh sorry, my bad for this one !!"

[Ajinkya]

Actually I am referring one textbook named "POWER SYSTEM DYNAMICS

AND STABILITY by Peter W. Sauer and M. A. Pai " in which they have modelled this IEEE 9bus system.

There they have taken [𝐼𝑑1 𝐼𝑞1 𝐼𝑑2 𝐼𝑞2 𝐼𝑑3 𝐼𝑞3] as state variables as well along with the

bus voltage magnitudes [V1 ... V9] and its angles [th1 ... th9].

Hence it turns out to be 45 state variables with 45 equations. (21 differential and 24 algebraic)

[Torsten]

Ok. Then modify the start of your function "fun_DAE" like

function dydt = fun_DAE(t,y,data1,data2,...)
    % Differential variables
    Eq1 = y(1);
    Ed1 = y(2);
    delta1 = y(3);
    omega1 = y(4);
    Efd1 = y(5);
    Rf1 = y(6);
    Vr1 = y(7);
    Eq2 = y(8);
    Ed2 = y(9);
    delta2 = y(10);
    omega2 = y(11);
    Efd2 = y(12);
    Rf2 = y(13);
    Vr2 = y(14); 
    Eq3 = y(15);
    Ed3 = y(16);
    delta3 = y(17);
    omega3 = y(18);
    Efd3 = y(19);
    Rf3 = y(20);
    Vr3 = y(21);  
    %Algebraic variables
    Id1 = y(22);
    Id2 = y(23);
    Id3 = y(24);
    Iq1 = y(25);
    Iq2 = y(26);
    Iq3 = y(27);
    V1 = y(28);
    V2 = y(29);
    V3 = y(30);
    V4 = y(31);
    V5 = y(32);
    V6 = y(33);
    V7 = y(34);
    V8 = y(35);
    V9 = y(36);
    theta1 = y(37);
    theta2 = y(38);
    theta3 = y(39);
    theta4 = y(40);
    theta5 = y(41);
    theta6 = y(42);
    theta7 = y(43);
    theta8 = y(44);
    theta9 = y(45);  
    ...
end

and program the right-hand sides of your differential and algebraic equations with the help of these variables Eq1,...,theta9.

Take care to supply the (45x1) vector of initial values in exactly this order.

Come back when you have modified "fun_DAE" as described above.

[Ajinkya]

Do I also need to provide the MassMatrix ?

[Torsten]

Yes. The mass matrix is

M = [eye(21),zeros(21,24);zeros(24,21),zeros(24,24)]

This assumes that the differential equations are all of the form

dy(i)/dt = rhs(i)

i.e. you have divided your equations by potential prefactors of dy(i)/dt (like T_do1, T_qo1, 2H1/omegas, T_E1, T_F1, T_A1, etc.).

[Ajinkya]

Warning: IDAS returned -4 from module IDAS function IDASolve: At t = 10.4731 and h = 8.61776e-09, the corrector convergence failed repeatedly or with |h| = hmin. 
S1 = 
  ODEResults with properties:

        Time: [0 9.7441e-05 1.9488e-04 3.8976e-04 7.7953e-04 0.0016 0.0031 0.0062 0.0094 0.0125 0.0187 0.0249 0.0312 0.0437 0.0561 0.0686 0.0811 0.1060 0.1310 0.1458 0.1607 … ] (1×431 double)
    Solution: [45×431 double]
>> 

When I use the idas solver then I am getting this warning.

[Torsten]

I don't know the code you are using at the moment. So the only thing that I can tell from the error message is that IDAS has a problem to integrate past t = 10.4731.

I suggest you plot and evaluate the results up to this time to see if the solution variables physically make sense.

[Ajinkya]

This is the code which when you run will give the warning

 clear all
    clc
    format short
    format compact 
% LOADS MAT FILE CONTAINING MACHINE DATA "dataM" and EXCITER DATA "dataE"
    load data.mat             
% THE ABOVE DATA IS DEFINED IN THE FILE "Machine_Exciter_Data.m"
% YBUS AND LOAD FLOW PROGRAM
    [Y,Vb] = YBus();
    I = Y*Vb;
    V = abs(Vb);     th = angle(Vb); 
% POWER INJECTIONS AT BUSES
    Sbus = Vb.*conj(I);
% CREATING A STRUCTURE NAMED 'BUS' TO STORE VALUES OF V, th, Y
    save Bus V Y th
%%
% MACHINE DATA
   [D_M, H, M, Xd, Xdd] = deal(data.D_M, data.H, data.M, data.Xd,data.Xdd);
   [Xq, Xqq, Td0, Tq0] = deal(data.Xq, data.Xqq, data.Td0, data.Tq0);
   D = D_M.*M;
   ws = 2*pi*60;
% EXCITER DATA
   [KA, TA, KE, TE, KF, TF] = deal(data.KA, data.TA, data.KE, data.TE, data.KF, data.TF);

% DETERMINING THE INITIAL VALUES OF THE STATE VARIABLES
for k = 1:3
    Pg(k) = real(Sbus(k));
    Qg(k) = imag(Sbus(k));
    Ig(k) = conj((Pg(k) + 1i*Qg(k))/(V(k)*exp(1i*th(k))));
    delt(k) = angle(V(k)*exp(1i*th(k)) + 1i*Xq(k)*Ig(k));
    Idq(k) = abs(Ig(k))*exp(1i*(angle(Ig(k))-delt(k) + pi/2));
    Id(k) = real(Idq(k));   Iq(k) = imag(Idq(k));
    Vdq(k) = V(k)*exp(1i*(th(k) - delt(k) + pi/2));
    Ed(k) = real(Vdq(k) - Xqq(k)*imag(Idq(k)));
    Eq(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k));
    Efd(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k));
    Rf(k) = KF(k)*(imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k)))/TF(k);
    Vr(k) = (KE(k) + 0.0039*exp(1.555*Efd(k)))*Efd(k);
    Vref(k) = V(k) + Vr(k)/KA(k);
    Tm(k) = Ed(k)*real(Idq(k)) + Eq(k)*imag(Idq(k)) + (Xqq(k) - Xdd(k))*real(Idq(k))*imag(Idq(k));
    w(k) = ws;
end
%%
for i = 1:3     % STATE VARIABLE VECTOR 
    %           1       2       3           4       5       6       7
    x(:,i) = [Eq(i); Ed(i); delt(i);    w(i);   Efd(i);    Rf(i); Vr(i)];
end
% InitVal STRUCT TO STORE THE INITIAL VALUES OF THE DIFFERENTIAL VARIABLE
    y0 = [x(:,1);x(:,2);x(:,3);Id(1,:)';Iq(1,:)';V(:,1);th(:,1)];
    Init.x = x;
    Init.Idq = [Id;Iq]; 
    Init.Vref_Tm = [Vref; Tm];
    Init.V = V;
    Init.th = th;
    Init.Y = Y;
    save y y0 
    save Init Init            % SAVING THE INITIAL VALUE STRUCT IN DATA FORMAT   
    MM = [eye(21),zeros(21,24);zeros(24,21),zeros(24,24)];
%%  SOLVING THE DAE
    F = ode;
    F.ODEFcn = @(t,y) func(t,y,data,Init);
    F.InitialValue = y0;
    F.MassMatrix = MM;
    F.Solver = 'idas';
    S1 = solve(F,0,100)

Warning: IDAS returned -4 from module IDAS function IDASolve: At t = 10.2132 and h = 1.2895e-08, the corrector convergence failed repeatedly or with |h| = hmin.

S1 = 
  ODEResults with properties:

        Time: [0 9.7441e-05 1.9488e-04 3.8976e-04 7.7953e-04 0.0016 0.0031 0.0062 0.0094 0.0125 0.0187 0.0249 0.0312 0.0437 0.0561 0.0686 0.0811 0.1060 0.1310 ... ] (1x366 double)
    Solution: [45x366 double]

    writematrix(S1.Time,'result.xlsx','Sheet','Time');
    writematrix(S1.Solution,'result.xlsx','Sheet','State')

[Torsten]

Replace

F.Solver = 'idas';

by

F.Solver = 'ode15s';

Further, it would be much more convenient for debugging and comparing your equations with those listed in your documentation if you would replace the y(:) in "func.m" by their variable names. That's why I rewrote the y-vector in local variables when I posted

    % Differential variables
    Eq1 = y(1);
    Ed1 = y(2);
    delta1 = y(3);
    omega1 = y(4);
    Efd1 = y(5);
    Rf1 = y(6);
    Vr1 = y(7);
    Eq2 = y(8);
    Ed2 = y(9);
    delta2 = y(10);
    omega2 = y(11);
    Efd2 = y(12);
    Rf2 = y(13);
    Vr2 = y(14); 
    Eq3 = y(15);
    Ed3 = y(16);
    delta3 = y(17);
    omega3 = y(18);
    Efd3 = y(19);
    Rf3 = y(20);
    Vr3 = y(21);  
    %Algebraic variables
    Id1 = y(22);
    Id2 = y(23);
    Id3 = y(24);
    Iq1 = y(25);
    Iq2 = y(26);
    Iq3 = y(27);
    V1 = y(28);
    V2 = y(29);
    V3 = y(30);
    V4 = y(31);
    V5 = y(32);
    V6 = y(33);
    V7 = y(34);
    V8 = y(35);
    V9 = y(36);
    theta1 = y(37);
    theta2 = y(38);
    theta3 = y(39);
    theta4 = y(40);
    theta5 = y(41);
    theta6 = y(42);
    theta7 = y(43);
    theta8 = y(44);
    theta9 = y(45);  
    

And you forgot to divide the right-hand sides of the differential equations by the multiplicative factors of the dy(i)/dt term (T_do1, T_qo1, 2H1/omegas, T_E1, T_F1, T_A1, etc.) (at least partially). E.g. Efd1 from the first equation was not divided by T_do1 - I didn't check what went wrong in the other equations.

[Ajinkya]

Okay, I will check it once and come back.

[Ajinkya]

Warning: Failure at t=3.927051e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(7.105427e-15) at time t. 
> In ode15s (line 625)
In matlab.ode.internal/LegacySolver/solve (line 140)
In ode/solveBetween (line 879)
In ode/solve (line 342)
In MainCode (line 64) 
S1 = 
  ODEResults with properties:

        Time: [0 3.3051e-05 6.6102e-05 9.9153e-05 4.2966e-04 7.6017e-04 0.0011 0.0044 0.0077 0.0110 0.0143 0.0245 0.0346 0.0448 … ] (1×104 double)
    Solution: [45×104 double]
>> 

again the similar warning, I have checked the equations again.

The solver is set to 'ode15s'.

[Torsten]

again the similar warning, I have checked the equations again.

Then you have an older MATLAB version or you changed something in the your equations what I told you was wrong. The code from above runs under MATLAB R2024b :

 clear all
    clc
    format short
    format compact 
% LOADS MAT FILE CONTAINING MACHINE DATA "dataM" and EXCITER DATA "dataE"
    load data.mat             
% THE ABOVE DATA IS DEFINED IN THE FILE "Machine_Exciter_Data.m"
% YBUS AND LOAD FLOW PROGRAM
    [Y,Vb] = YBus();
    I = Y*Vb;
    V = abs(Vb);     th = angle(Vb); 
% POWER INJECTIONS AT BUSES
    Sbus = Vb.*conj(I);
% CREATING A STRUCTURE NAMED 'BUS' TO STORE VALUES OF V, th, Y
    save Bus V Y th
%%
% MACHINE DATA
   [D_M, H, M, Xd, Xdd] = deal(data.D_M, data.H, data.M, data.Xd,data.Xdd);
   [Xq, Xqq, Td0, Tq0] = deal(data.Xq, data.Xqq, data.Td0, data.Tq0);
   D = D_M.*M;
   ws = 2*pi*60;
% EXCITER DATA
   [KA, TA, KE, TE, KF, TF] = deal(data.KA, data.TA, data.KE, data.TE, data.KF, data.TF);
% DETERMINING THE INITIAL VALUES OF THE STATE VARIABLES
for k = 1:3
    Pg(k) = real(Sbus(k));
    Qg(k) = imag(Sbus(k));
    Ig(k) = conj((Pg(k) + 1i*Qg(k))/(V(k)*exp(1i*th(k))));
    delt(k) = angle(V(k)*exp(1i*th(k)) + 1i*Xq(k)*Ig(k));
    Idq(k) = abs(Ig(k))*exp(1i*(angle(Ig(k))-delt(k) + pi/2));
    Id(k) = real(Idq(k));   Iq(k) = imag(Idq(k));
    Vdq(k) = V(k)*exp(1i*(th(k) - delt(k) + pi/2));
    Ed(k) = real(Vdq(k) - Xqq(k)*imag(Idq(k)));
    Eq(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k));
    Efd(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k));
    Rf(k) = KF(k)*(imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k)))/TF(k);
    Vr(k) = (KE(k) + 0.0039*exp(1.555*Efd(k)))*Efd(k);
    Vref(k) = V(k) + Vr(k)/KA(k);
    Tm(k) = Ed(k)*real(Idq(k)) + Eq(k)*imag(Idq(k)) + (Xqq(k) - Xdd(k))*real(Idq(k))*imag(Idq(k));
    w(k) = ws;
end
%%
for i = 1:3     % STATE VARIABLE VECTOR 
    %           1       2       3           4       5       6       7
    x(:,i) = [Eq(i); Ed(i); delt(i);    w(i);   Efd(i);    Rf(i); Vr(i)];
end
% InitVal STRUCT TO STORE THE INITIAL VALUES OF THE DIFFERENTIAL VARIABLE
    y0 = [x(:,1);x(:,2);x(:,3);Id(1,:)';Iq(1,:)';V(:,1);th(:,1)];
    Init.x = x;
    Init.Idq = [Id;Iq]; 
    Init.Vref_Tm = [Vref; Tm];
    Init.V = V;
    Init.th = th;
    Init.Y = Y;
    save y y0 
    save Init Init            % SAVING THE INITIAL VALUE STRUCT IN DATA FORMAT   
    MM = [eye(21),zeros(21,24);zeros(24,21),zeros(24,24)];
%%  SOLVING THE DAE
    F = ode;
    F.ODEFcn = @(t,y) func(t,y,data,Init);
    F.InitialValue = y0;
    F.MassMatrix = MM;
    F.Solver = 'ode15s';
    S1 = solve(F,0,100)

S1 = 
  ODEResults with properties:

        Time: [0 8.9788e-05 1.7958e-04 2.6937e-04 0.0012 0.0021 0.0030 0.0105 0.0181 0.0257 0.0333 0.0514 0.0695 0.0876 0.1057 0.1238 0.1436 0.1634 0.1833 ... ] (1x3762 double)
    Solution: [45x3762 double]

[Ajinkya]

I have the 2024b version installed.

In the MassMatrix, do I need to mention that it is singular. I think there is an option to do that.

Among these properties, can I set the Jacobian matrix? I already have the matrix with me.

>> F
F = 
  ode with properties:

   Problem definition
               ODEFcn: @(t,y)func(t,y,data,Init)
          InitialTime: 0
         InitialValue: [45×1 double]
             Jacobian: []
           MassMatrix: [1×1 odeMassMatrix]
         EquationType: standard
   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 5.0000e-04
               Solver: ode15s
        SolverOptions: [1×1 matlab.ode.options.ODE15s]
  Show all properties

            EquationType: standard
                  ODEFcn: @(t,y)func(t,y,data,Init)
              Parameters: []
             InitialTime: 0
            InitialValue: [45×1 double]
            InitialSlope: [0×1 double]
                Jacobian: []
              MassMatrix: [1×1 odeMassMatrix]
    NonNegativeVariables: []
         EventDefinition: []
             Sensitivity: []
       AbsoluteTolerance: 1.0000e-06
       RelativeTolerance: 5.0000e-04
                  Solver: ode15s
           SolverOptions: [1×1 matlab.ode.options.ODE15s]
          SelectedSolver: ode15s
>> 

[Torsten]

I have the 2024b version installed.

What then was the reason for the error message you got ? Did you make changes to the code ?

In the MassMatrix, do I need to mention that it is singular. I think there is an option to do that.

Just try if it changes anything. I don't think so.

Among these properties, can I set the Jacobian matrix? I already have the matrix with me.

Are you sure ? You computed the derivatives of the expressions you defined in "fun_DAE" with respect to all the y-variables ? I strongly doubt this is true. In case you really have this monster: Specify the Jacobian property as a handle to a function that computes the matrix elements and that accepts two input arguments, dfdy = Fjac(t,y).

[Ajinkya]

I am doing some changes in that function of Jacobian. I have defined it as

function Jcb = jacob(y,data)

end

The Jacobian will turn up to be of 45x45. Is the function handle is correct or should I make some changes.

[Torsten]

The body of the function is correct. But I don't think it's worth spending your time with: even if you don't make errors in setting up this function, it will not give you better results. Only the computation time might be reduced by a millisecond.

[Ajinkya]

Okay. I have added the code in the below "More Answers" section.

I am getting the results but they are not satisfactory

[Torsten]

Implementing such a complicated DAE system does not mean writing the equations, clicking on the RUN button and getting results as expected. It's an iterative process: you will have to check the results, make a guess what could be wrong in the code when you look at the curves, check and perhaps modify parameters and equations again, run the code, ... .

I did what I could to make the code work technically. Now it's up to you to make it give results that physically make sense.

[Ajinkya]

Surely. I will do the needful.

Thank you so much for your valuable suggestions and help.

I was stuck with this problem for days and now I have done something fruitful.

Grateful for the time you have given !!

Answer 1

[Steven Lord]

Since you're using release R2024b, you have access to three of the SUNDIALS solvers via the ode object. This functionality was added in release R2024a as stated in the Release Notes here. Try using the object to set up your problem and solve it using "idas" as the value of the Solver property.

[Ajinkya]

could you please guide me how do I proceed in this?

[Torsten]

F = ode;
F.ODEFcn = @(t,y) 2*t;
F.InitialValue = 0;
F.Solver = 'idas';
F

F = 
  ode with properties:

   Problem definition
               ODEFcn: @(t,y)2*t
          InitialTime: 0
         InitialValue: 0
             Jacobian: []
         EquationType: standard

   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 1.0000e-03
               Solver: idas

  Show all properties

sol = solve(F,0,10);
plot(sol.Time,sol.Solution)

●

Answer 2

[Ajinkya]

    clear 
    clc
    format short
    format compact 
% LOADS MAT FILE CONTAINING MACHINE DATA "dataM" and EXCITER DATA "dataE"
    load data        
% THE ABOVE DATA IS DEFINED IN THE FILE "Machine_Exciter_Data.m"
% YBUS AND LOAD FLOW PROGRAM
    [Y,Vb] = YBus();
    I = Y*Vb;
    V = abs(Vb);     th = angle(Vb); 
    % POWER INJECTIONS AT BUSES
    Sbus = Vb.*conj(I);
    % CREATING A STRUCTURE NAMED 'BUS' TO STORE VALUES OF V, th, Y
    save Bus V Y th
%%
    % MACHINE DATA
   [D_M, H, M, Xd, Xdd] = deal(data.D_M, data.H, data.M, data.Xd,data.Xdd);
   [Xq, Xqq, Td0, Tq0] = deal(data.Xq, data.Xqq, data.Td0, data.Tq0);
   D = D_M.*M;
   ws = 2*pi*60;
    % EXCITER DATA
   [KA, TA, KE, TE, KF, TF] = deal(data.KA, data.TA, data.KE, data.TE, data.KF, data.TF);

   % DETERMINING THE INITIAL VALUES OF THE STATE VARIABLES
   for k = 1:3
    Pg(k) = real(Sbus(k));
    Qg(k) = imag(Sbus(k));
    Ig(k) = conj((Pg(k) + 1i*Qg(k))/(V(k)*exp(1i*th(k))));
    delt(k) = angle(V(k)*exp(1i*th(k)) + 1i*Xq(k)*Ig(k));
    Idq(k) = abs(Ig(k))*exp(1i*(angle(Ig(k))-delt(k) + pi/2));
    Id(k) = real(Idq(k));   Iq(k) = imag(Idq(k));
    Vdq(k) = V(k)*exp(1i*(th(k) - delt(k) + pi/2));
    Ed(k) = real(Vdq(k) - Xqq(k)*imag(Idq(k)));
    Eq(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k));
    Efd(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k));
    Rf(k) = KF(k)*(imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k)))/TF(k);
    Vr(k) = (KE(k) + 0.0039*exp(1.555*Efd(k)))*Efd(k);
    Vref(k) = V(k) + Vr(k)/KA(k);
    Tm(k) = Ed(k)*real(Idq(k)) + Eq(k)*imag(Idq(k)) + (Xqq(k) - Xdd(k))*real(Idq(k))*imag(Idq(k));
    w(k) = ws;
    end
%%
    for i = 1:3     % STATE VARIABLE VECTOR 
        %           1       2       3           4       5       6       7
        x(:,i) = [Eq(i); Ed(i); delt(i);    w(i);   Efd(i);    Rf(i); Vr(i)];
    end
    % InitVal STRUCT TO STORE THE INITIAL VALUES OF THE DIFFERENTIAL VARIABLE
    y0 = [x(:,1);x(:,2);x(:,3);Id(1,:)';Iq(1,:)';V(:,1);th(:,1)];
    Init.x = x;
    Init.Idq = [Id;Iq]; 
    Init.Vref = Vref;   Init.Tm = Tm;
    Init.V = V;
    Init.th = th;
    Init.Y = Y;
    save y y0 
    save Init Init            % SAVING THE INITIAL VALUE STRUCT IN DATA FORMAT   
    
    % Eq1 = y(1);             Eq2 = y(8);             Eq3 = y(15);
    % Ed1 = y(2);             Ed2 = y(9);             Ed3 = y(16);
    % delta1 = y(3);          delta2 = y(10);         delta3 = y(17);
    % omega1 = y(4);          omega2 = y(11);         omega3 = y(18);
    % Efd1 = y(5);            Efd2 = y(12);           Efd3 = y(19);
    % Rf1 = y(6);             Rf2 = y(13);            Rf3 = y(20);
    % Vr1 = y(7);             Vr2 = y(14);            Vr3 = y(21);
    % 
    %  Id1 = y(22);           V1 = y(28);             theta1 = y(37);
    %  Id2 = y(23);           V2 = y(29);             theta2 = y(38);
    %  Id3 = y(24);           V3 = y(30);             theta3 = y(39);
    %  Iq1 = y(25);           V4 = y(31);             theta4 = y(40);
    %  Iq2 = y(26);           V5 = y(32);             theta5 = y(41);
    %  Iq3 = y(27);           V6 = y(33);             theta6 = y(42);
    %                         V7 = y(34);             theta7 = y(43);
    %                         V8 = y(35);             theta8 = y(44);
    %                         V9 = y(36);             theta9 = y(45);
%%  SOLVING THE DAE
    F = ode;
    F.ODEFcn = @(t,y) func(t,y,data,Init);
    F.InitialValue = y0;
    F.MassMatrix = massMatrix();
    F.RelativeTolerance = 1e-03;
    F.Solver = 'ode15s';
    %F.Jacobian = @(t,y,data) jcb;
   % F.SolverOptions.MaxOrder = 3;   % NOT FOR ODE23x
    F.SolverOptions.MaxStep = 0.1;
    S1 = solve(F,0,100)
    % writematrix(S1.Time,'result.xlsx','Sheet','Time');
    % writematrix(S1.Solution,'result.xlsx','Sheet','State')
%%
    plot(S1.Time, S1.Solution(1, :), '-r', 'DisplayName', 'Eq1');
    hold on;
    plot(S1.Time, S1.Solution(8, :), '-g', 'DisplayName', 'Eq2');
    plot(S1.Time, S1.Solution(15, :), '-b', 'DisplayName', 'Eq3');
    hold off;
    xlabel('Time');
    ylabel('Flux Linkages');
    legend show;
    grid on;