Solve DAEs using IDASolve function of SUNDIALS 2.6.2 version

3 次查看（过去 30 天）
Ajinkya
Ajinkya 2025-4-22
评论： Ajinkya 2025-5-1
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.
  34 个评论
显示 32更早的评论
Torsten
Torsten 2025-5-1
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
Ajinkya 2025-5-1
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 !!

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Steven Lord
Steven Lord 2025-4-22
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.
  2 个评论
Torsten
Torsten 2025-4-23
Ran in:
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)
●

请先登录，再进行评论。

更多回答（1 个）

Ajinkya
Ajinkya 2025-5-1
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;

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by