Solve DAEs using IDASolve function of SUNDIALS 2.6.2 version
Show older comments
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.
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.
34 Comments
Torsten
on 23 Apr 2025
Could you include your code as an .m-file so that we are able to run it here directly ?
Torsten
on 26 Apr 2025
"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
on 26 Apr 2025
Torsten
on 26 Apr 2025
What are the unknowns in the algebraic equations ?
Ajinkya
on 27 Apr 2025
Torsten
on 27 Apr 2025
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.
Torsten
on 27 Apr 2025
"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
on 27 Apr 2025
Ajinkya
on 28 Apr 2025
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
on 28 Apr 2025
Torsten
on 28 Apr 2025
So there are additional variables in the algebraic equations that are unknown ? Or how to interprete "Oh sorry, my bad for this one !!"
Ajinkya
on 28 Apr 2025
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
on 28 Apr 2025
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
on 29 Apr 2025
Torsten
on 29 Apr 2025
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.
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
on 29 Apr 2025
Ajinkya
on 29 Apr 2025
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)
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).
Torsten
on 1 May 2025
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
on 1 May 2025
Torsten
on 1 May 2025
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
on 1 May 2025
Accepted Answer
Steven Lord
on 22 Apr 2025
0 votes
2 Comments
Ajinkya
on 23 Apr 2025
F = ode;
F.ODEFcn = @(t,y) 2*t;
F.InitialValue = 0;
F.Solver = 'idas';
F
sol = solve(F,0,10);
plot(sol.Time,sol.Solution)
More Answers (1)
Ajinkya
on 1 May 2025
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
